1D Transient Thermal Conduction, BCs of 1st Kind¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
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.6, pg.19, “A Transient 1D Temperature Distribution, Time-Dependent Boundary Conditions of 1st Kind.”
The domain is a 20x1x1 meter rectangular column extending along the positive and negative x-axis and is made up of 20x1x1 cubic grid cells with dimensions 1x1x1 meters. The domain material is assigned the following properties: thermal conductivity K = 2.0 W/(m-C); specific heat capacity Cp = 1.5 J/(m-C); density rho = 2,500 kg/m^3.
The temperature is initially uniform at T(t=0) = 0 C. The boundary temperatures are:
\[ \begin{align}\begin{aligned}T(-L,t) = T_b t\\T(L,t) = T_b t\end{aligned}\end{align} \]
where L = 10 m and \(T_b\) = 2 C/day. The transient temperature distribution is governed by,
\[\rho c_p {{\partial T} \over {\partial t}} = K {{\partial^{2} T} \over {\partial x^{2}}}\]
With the given boundary conditions, the solution is defined by,
\[ \begin{align}\begin{aligned}\chi = {K \over {\rho c_p}}\\T(x,t) = T_b t + {{T_b(x^2-L^2)}\over{2\chi}} + {{16T_bL^2}\over{\chi\pi^3}} \sum_{n=0}^{\infty}{{(-1)^n}\over{(2n+1)^3}} cos{\left({{(2n+1)x\pi}\over{2L}}\right)} e^{\left({-\chi(2n+1)^2\pi^2{{t}\over{4L^2}}}\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 20 1 1
DXYZ
1.0d0
1.0d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
20.0d0 1.0d0 1.0d0
/
END
REGION left_end
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END
REGION right_end
FACE EAST
COORDINATES
20.0d0 0.d0 0.d0
20.0d0 1.0d0 1.0d0
/
END
MATERIAL_PROPERTY beam
ID 1
POROSITY 1.d-20
TORTUOSITY 1.d0
ROCK_DENSITY 2.5d3 kg/m^3
HEAT_CAPACITY 1.5 J/kg-C
THERMAL_CONDUCTIVITY_DRY 2.0 W/m-C
THERMAL_CONDUCTIVITY_WET 2.0 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-13
PERM_Y 1.d-13
PERM_Z 1.d-13
/
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
EOS WATER
DENSITY CONSTANT 1000. kg/m^3
#ENTHALPY CONSTANT 0.d0 J/kmol
VISCOSITY CONSTANT 1.d-3 Pa-s
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
STRATA
REGION all
MATERIAL beam
END
TIME
FINAL_TIME 0.75 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
TIMES day 0.25 0.50
#NO_PRINT_INITIAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 0.0d0 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 LIST
# T = Tb*t; Tb=2C
TIME_UNITS day
DATA_UNITS C
INTERPOLATION LINEAR
#time #temperature
0.00d0 0.0d0
0.25d0 0.5d0
0.50d0 1.0d0
1.00d0 2.0d0
/
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION right_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE LIST
# T = Tb*t; Tb=2C
TIME_UNITS day
DATA_UNITS C
INTERPOLATION LINEAR
#time #temperature
0.00d0 0.0d0
0.25d0 0.5d0
0.50d0 1.0d0
1.00d0 2.0d0
/
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 20 1 1
DXYZ
1.0d0
1.0d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
20.0d0 1.0d0 1.0d0
/
END
REGION left_end
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END
REGION right_end
FACE EAST
COORDINATES
20.0d0 0.d0 0.d0
20.0d0 1.0d0 1.0d0
/
END
MATERIAL_PROPERTY beam
ID 1
POROSITY 1.d-20
TORTUOSITY 1.d0
ROCK_DENSITY 2.5d3 kg/m^3
SPECIFIC_HEAT 1.5 J/kg-C
THERMAL_CONDUCTIVITY_DRY 2.0 W/m-C
THERMAL_CONDUCTIVITY_WET 2.0 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-13
PERM_Y 1.d-13
PERM_Z 1.d-13
/
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 0.75 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
TIMES day 0.25 0.50
#NO_PRINT_INITIAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 0.0D0 C
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION left_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE LIST
# T = Tb*t; Tb=2 C/day
TIME_UNITS day
DATA_UNITS C
INTERPOLATION LINEAR
#time #temperature
0.00d0 0.00d0
0.25d0 0.50d0
0.50d0 1.00d0
1.00d0 2.00d0
/
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION right_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE LIST
# T = Tb*t; Tb=2C
TIME_UNITS day
DATA_UNITS C
INTERPOLATION LINEAR
#time #temperature
0.00d0 0.00d0
0.25d0 0.50d0
0.50d0 1.00d0
1.00d0 2.00d0
/
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
ONLY_ENERGY_EQ
END_SUBSURFACE
The Python Script¶
def thermal_transient_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.6, pg.19
# "A Transient 1D Temperature Distribution, Time-Dependent Boundary Conditions
# of 1st Kind"
# With some modification of the parameter values given in Kolditz (2015)
#
# Author: Jenn Frederick
# Date: 06/30/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] = 20. # [m] lx
nxyz[0] = 10 # [m] nx
dxyz = lxyz/nxyz # [m]
K = 2.0 # [W/m-C]
Cp = 1.5 # [J/kg-C]
rho = 2500. # [kg/m^3]
Tb_day = 2.0 # [C/day]
Tb_sec = 2.0/(24.0*3600.0) # [C/sec]
chi = K/(Cp*rho)
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_soln1 = np.linspace(0.+(dx/2.),(Lx/2)-(dx/2.),int(nx/2)) # [m]
x_soln2 = np.linspace(-((Lx/2)-(dx/2.)),-(0.+(dx/2.)),int(nx/2)) # [m]
x_soln = np.concatenate((x_soln2,x_soln1),axis=0) # [m]
# Add boundary values to analytical solution
x_soln = np.concatenate(([-Lx/2],x_soln),axis=0) # [m]
x_soln = np.concatenate((x_soln,[Lx/2]),axis=0) # [m]
t_soln = np.array([0.0,0.25,0.50,0.75]) # [day]
T_soln = np.zeros((4,(nx+2))) # [C]
x_pflotran1 = np.linspace(0.+(dx/2.),(Lx/2)-(dx/2.),int(nx/2)) # [m]
x_pflotran2 = np.linspace(-((Lx/2)-(dx/2.)),-(0.+(dx/2.)),int(nx/2)) # [m]
x_pflotran = np.concatenate((x_pflotran2,x_pflotran1),axis=0) # [m]
T_pflotran = np.zeros((4,nx)) # [C]
# create the analytical solution
for time in range(4):
t = t_soln[time]*(24.0*3600.0) # [sec]
T_soln[time,:] = Tb_sec*t + ((Tb_sec*(pow(x_soln,2)-pow(Lx/2,2)))/(2.*chi))
sum_term = np.zeros(nx+2)
sum_term_old = np.zeros(nx+2)
n = 0
epsilon = 1.0
while epsilon > epsilon_value:
sum_term_old = sum_term
sum_term = sum_term_old + (((pow(-1.,n))/(pow(((2*n)+1),3)))*np.cos((math.pi*x_soln*((2*n)+1))/(2*Lx/2))*np.exp(-chi*pow((2*n)+1,2)*pow(math.pi,2)*(t/(4*pow(Lx/2,2)))))
epsilon = np.max(np.abs(sum_term_old-sum_term))
n = n + 1
T_soln[time,:] = T_soln[time,:] + ((16.*Tb_sec*pow(Lx/2,2))/(chi*pow(math.pi,3)))*sum_term
# run PFLOTRAN simulation
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)
# load data from HDF5
index_string = 'Time: 0.00000E+00 d/Temperature [C]'
T_pflotran[0,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_1st_kind.h5',index_string,False)
index_string = 'Time: 2.50000E-01 d/Temperature [C]'
T_pflotran[1,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_1st_kind.h5',index_string,False)
index_string = 'Time: 5.00000E-01 d/Temperature [C]'
T_pflotran[2,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_1st_kind.h5',index_string,False)
index_string = 'Time: 7.50000E-01 d/Temperature [C]'
T_pflotran[3,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_1st_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_transient(path,t_soln,'day',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 Transient Thermal Conduction, BCs of 2nd Kind¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
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.7, pg.20, “A Transient 1D Temperature Distribution, Time-Dependent Boundary Conditions of 2nd Kind.”
The domain is a 25x1x1 meter rectangular column extending along the positive x-axis and is made up of 50x1x1 cubic grid cells with dimensions 0.5x1x1 meters. The domain material is assigned the following properties: thermal conductivity K = 1.16 W/(m-C); specific heat capacity Cp = 0.01 J/(m-C); density \(\rho\) = 2,000 kg/m^3. The accuracy of the numerical solution is very sensitive to the time step size because the heat flux value is held constant at the value at the beginning of the timestep. Reducing the timestep will make the numerical solution more accurate relative to the analytical solution, but note this is not an issue with time truncation error.
The temperature is initially uniform at T(t=0) = 0 C. The heat flux boundary conditions (in units of W/m2) are:
\[ \begin{align}\begin{aligned}q(0,t) = 0\\q(L,t) = 0.385802 t\end{aligned}\end{align} \]
where L = 25 m. The transient temperature distribution is governed by,
\[\rho c_p {{\partial T} \over {\partial t}} = K {{\partial^{2} T} \over {\partial x^{2}}}\]
With the given boundary conditions, the solution is defined by,
\[T(x,t) = {{8q\sqrt{\chi t^3}} \over K} \sum_{n=0}^{\infty} \left[{ i^3erfc{{(2n+1)L-x}\over{2\sqrt{\chi t}}} + i^3erfc{{(2n+1)L+x}\over{2\sqrt{\chi t}}} }\right]\]
where \(i^3erfc(g)\) represents the third repeated integral of the complimentary error function, given by,
\[i^3erfc(g) = {2 \over \pi} \int^{\infty}_{g} {{(s-g)^3}\over{3!}} e^{-s^2} ds\]
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
2.5d0
1.0d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
25.0d0 1.0d0 1.0d0
/
END
REGION left_end
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END
REGION right_end
FACE EAST
COORDINATES
25.0d0 0.d0 0.d0
25.0d0 1.0d0 1.0d0
/
END
MATERIAL_PROPERTY beam
ID 1
POROSITY 1.d-10
TORTUOSITY 1.d0
ROCK_DENSITY 2.0d3 kg/m^3
HEAT_CAPACITY 0.01 J/kg-C
THERMAL_CONDUCTIVITY_DRY 1.16 W/m-C
THERMAL_CONDUCTIVITY_WET 1.16 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-15
PERM_Y 1.d-15
PERM_Z 1.d-15
/
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
EOS WATER
DENSITY CONSTANT 1000. kg/m^3
#ENTHALPY CONSTANT 1.d0 J/kmol
VISCOSITY CONSTANT 1.d-3 Pa-s
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
NUMERICAL_METHODS FLOW
LINEAR_SOLVER
SOLVER DIRECT
LU_ZERO_PIVOT_TOL 1.d-20
/
END
STRATA
REGION all
MATERIAL beam
END
TIME
FINAL_TIME 0.20 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
TIMES day 0.01 0.04 0.09 0.12
NO_PRINT_INITIAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 0.d0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION left_end
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX 0.d0 W/m^2
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END
FLOW_CONDITION right_end
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX LIST
# flux = 0.385802*t [W/m^2]
TIME_UNITS day
DATA_UNITS W/m^2
INTERPOLATION LINEAR
#time #heat_flux
0 0
0.12 0.04629624
/
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
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 10 1 1
DXYZ
2.5d0
1.0d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
25.0d0 1.0d0 1.0d0
/
END
REGION left_end
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END
REGION right_end
FACE EAST
COORDINATES
25.0d0 0.d0 0.d0
25.0d0 1.0d0 1.0d0
/
END
MATERIAL_PROPERTY beam
ID 1
POROSITY 1.d-10
TORTUOSITY 1.d0
ROCK_DENSITY 2.0d3 kg/m^3
SPECIFIC_HEAT 0.01 J/kg-C
THERMAL_CONDUCTIVITY_DRY 1.16 W/m-C
THERMAL_CONDUCTIVITY_WET 1.16 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
EOS WATER
DENSITY CONSTANT 1000. kg/m^3
#ENTHALPY CONSTANT 1.d0 J/kmol
VISCOSITY CONSTANT 1.d-3 Pa-s
END
NUMERICAL_METHODS FLOW
LINEAR_SOLVER
SOLVER DIRECT
LU_ZERO_PIVOT_TOL 1.d-20
/
END
STRATA
REGION all
MATERIAL beam
END
TIME
FINAL_TIME 0.20 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.0001 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
TIMES day 0.01 0.04 0.09 0.12
NO_PRINT_INITIAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 0.0d0 C
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION left_end
TYPE
ENERGY_FLUX neumann
LIQUID_PRESSURE dirichlet
/
ENERGY_FLUX 0.d0 W/m^2
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION right_end
TYPE
ENERGY_FLUX neumann
LIQUID_PRESSURE dirichlet
/
ENERGY_FLUX LIST
# flux = 0.385802*t [W/m^2]
TIME_UNITS day
DATA_UNITS W/m^2
INTERPOLATION LINEAR
#time #heat_flux
0.00d0 0.0d0
1.00d0 0.385802
/
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
ONLY_ENERGY_EQ
END_SUBSURFACE
The Python Script¶
def thermal_transient_1D_BC2ndkind(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.7, pg.20
# "A Transient 1D Temperature Distribution, Time-Dependent Boundary Conditions
# of 2nd Kind"
# With some modification of the parameter values given in Kolditz (2015)
#
# Author: Jenn Frederick
# Date: 07/12/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] = 25. # [m] lx
nxyz[0] = 10 # [m] nx
dxyz = lxyz/nxyz # [m]
K = 1.16 # [W/m-C]
Cp = 0.01 # [J/kg-C]
rho = 2000. # [kg/m^3]
q = 0.385802 # [C/day]
chi = K/(Cp*rho)
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_soln = np.concatenate((x_soln,[Lx]),axis=0) # [m]
T_soln = np.zeros((4,nx+1)) # [C]
x_pflotran = x_soln[0:nx] # [m]
T_pflotran = np.zeros((4,nx)) # [C]
t_soln = np.array([0.01,0.04,0.09,0.12]) # [day]
# create the analytical solution
p1 = np.zeros(nx+1)
p2 = np.zeros(nx+1)
i3erfc1 = np.zeros(nx+1)
i3erfc2 = np.zeros(nx+1)
for time in range(4):
t = t_soln[time]*24.0*3600.0 # [sec]
sum_term = np.zeros(nx+1)
sum_term_old = np.zeros(nx+1)
n = 0
epsilon = 1.0
while epsilon > epsilon_value:
p1 = ((((2*n)+1)*Lx)-x_soln)/(2.*math.sqrt(chi*t))
i3erfc1 = (1./12.)*np.sqrt(math.pi)*p1**3*erf(1.0*p1) - (1./12.)*np.sqrt(math.pi)*p1**3 + (1./3.)*p1**2*np.exp(-1.0*p1**2) + 0.5*p1*(-0.5*p1*np.exp(-1.0*p1**2) + 0.25*np.sqrt(math.pi)*erf(1.0*p1)) - 0.125*np.sqrt(math.pi)*p1 + (1./12.)*np.exp(-1.0*p1**2)
p2 = ((((2*n)+1)*Lx)+x_soln)/(2.*math.sqrt(chi*t))
i3erfc2 = (1./12.)*np.sqrt(math.pi)*p2**3*erf(1.0*p2) - (1./12.)*np.sqrt(math.pi)*p2**3 + (1./3.)*p2**2*np.exp(-1.0*p2**2) + 0.5*p2*(-0.5*p2*np.exp(-1.0*p2**2) + 0.25*np.sqrt(math.pi)*erf(1.0*p2)) - 0.125*np.sqrt(math.pi)*p2 + (1./12.)*np.exp(-1.0*p2**2)
sum_term_old = sum_term
sum_term = sum_term_old + (2./math.sqrt(math.pi))*i3erfc1 + (2./math.sqrt(math.pi))*i3erfc2
epsilon = np.max(np.abs(sum_term_old-sum_term))
n = n + 1
T_soln[time,:] = ((8*q*(1./(24.*3600.))*math.sqrt(chi*pow(t,3)))/K)*sum_term
# To calculate the definite integral, use sympy:
# p1, s = sym.symbols('p1 s')
# i3erfc1 = sym.integrate(((pow((s-p1),3))/(3.*2.*1.))*sym.exp(-1.*pow(s,2)),(s,p1,sym.oo))
# p2, s = sym.symbols('p2 s')
# i3erfc2 = sym.integrate(((pow((s-p2),3))/(3.*2.*1.))*sym.exp(-1.*pow(s,2)),(s,p2,sym.oo))
# run PFLOTRAN simulation
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)
index_string = 'Time: 1.00000E-02 d/Temperature [C]'
T_pflotran[0,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_2nd_kind.h5',index_string,False)
index_string = 'Time: 4.00000E-02 d/Temperature [C]'
T_pflotran[1,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_2nd_kind.h5',index_string,False)
index_string = 'Time: 9.00000E-02 d/Temperature [C]'
T_pflotran[2,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_2nd_kind.h5',index_string,False)
index_string = 'Time: 1.20000E-01 d/Temperature [C]'
T_pflotran[3,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_2nd_kind.h5',index_string,remove)
ierr = check(T_pflotran)
max_percent_error = calc_relative_error(T_soln[:,0:nx],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_transient(path,t_soln,'day',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 Transient Thermal Conduction, BCs of 1st & 2nd Kind¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
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.8, pg.21, “A Transient 1D Temperature Distribution, Non-Zero Initial Temperature, Boundary Conditions of 1st and 2nd Kind.”
The domain is a 100x1x1 meter rectangular beam extending along the positive x-axis and is made up of 50x1x1 hexahedral grid cells with dimensions 2x1x1 meters. The domain is composed of a single material and is assigned the following properties: thermal conductivity K = 0.5787037 W/(m-C); specific heat capacity Cp = 0.01 J/(m-C); density \(\rho\) = 2,000 kg/m^3.
The temperature is initially distributed according to T(x,t=0)=f(x), where f(x) is defined as
\[ \begin{align}\begin{aligned}f(x) = 0 \hspace{0.25in} 0 \leq x < {L \over 10}\\f(x) = {{10x} \over {3L}}-{1 \over 3} \hspace{0.25in} {L \over 10} \leq x < {{4L} \over 10}\\f(x) = 1 \hspace{0.25in} {{4L} \over 10} \leq x < {{6L} \over 10}\\f(x) = 3-{{10x} \over {3L}} \hspace{0.25in} {{6L} \over 10} \leq x < {{9L} \over 10}\\f(x) = 0 \hspace{0.25in} {{9L} \over 10} \leq x \leq L\end{aligned}\end{align} \]
At the two boundaries, a no heat flux condition is applied,
\[ \begin{align}\begin{aligned}q(0,t) = 0\\q(L,t) = 0\end{aligned}\end{align} \]
where L = 100 m. The transient temperature distribution is governed by,
\[\rho c_p {{\partial T} \over {\partial t}} = K {{\partial^{2} T} \over {\partial x^{2}}}\]
With the initial temperature given, the solution is defined by,
\[ \begin{align}\begin{aligned}T(x,t) = {1 \over 2} + \sum_{n=1}^{\infty} exp\left({-\chi n^2 \pi^2 {t \over L^2}}\right)\left({80 \over {3(n\pi)^2}}\right) cos{{n \pi y} \over L} cos{{n\pi} \over 2} sin{{n\pi} \over 4} sin{{3n\pi} \over 20}\\\chi = {K \over {\rho c_p}}\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 40 1 1
DXYZ
2.5d0
1.0d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 1.0d0 1.0d0
/
END
REGION left_end
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END
REGION right_end
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 1.0d0 1.0d0
/
END
MATERIAL_PROPERTY beam
ID 1
POROSITY 1.d-15
TORTUOSITY 1.d0
ROCK_DENSITY 2.0d3 kg/m^3
HEAT_CAPACITY 0.01 J/kg-C
THERMAL_CONDUCTIVITY_DRY 0.5787037 W/m-C
THERMAL_CONDUCTIVITY_WET 0.5787037 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-25
PERM_Y 1.d-25
PERM_Z 1.d-25
/
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 0.50 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.001 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
TIMES day 0.015 0.05 0.1
#NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET temperature_initial
HDF5_DATASET_NAME initial
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_initial
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION left_end
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX 0.d0 W/m^2
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END
FLOW_CONDITION right_end
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX 0.d0 W/m^2
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
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 40 1 1
DXYZ
2.5d0
1.0d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 1.0d0 1.0d0
/
END
REGION left_end
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END
REGION right_end
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 1.0d0 1.0d0
/
END
MATERIAL_PROPERTY beam
ID 1
POROSITY 1.d-10
TORTUOSITY 1.d0
ROCK_DENSITY 2.0d3 kg/m^3
SPECIFIC_HEAT 0.01 J/kg-C
THERMAL_CONDUCTIVITY_DRY 0.5787037 W/m-C
THERMAL_CONDUCTIVITY_WET 0.5787037 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 0.50 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.001 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
TIMES day 0.015 0.05 0.1
#NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET temperature_initial
HDF5_DATASET_NAME initial
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_initial
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION left_end
TYPE
ENERGY_FLUX neumann
LIQUID_PRESSURE dirichlet
/
ENERGY_FLUX 0.d0 W/m^2
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION right_end
TYPE
ENERGY_FLUX neumann
LIQUID_PRESSURE dirichlet
/
ENERGY_FLUX 0.d0 W/m^2
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
ONLY_ENERGY_EQ
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 initial condition
L = 100.
h5grp = h5file.create_group('initial')
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 = 101
rarray = np.zeros(nx,'=f8')
for i in range(nx):
if (0. <= i < (L/10.)):
rarray[i] = 0.
if ((L/10.) <= i < (4.*L/10.)):
rarray[i] = (10./(3.*L))*float(i) - (1./3.)
if ((4.*L/10.) <= i < (6.*L/10.)):
rarray[i] = 1.
if ((6.*L/10.) <= i < (9.*L/10.)):
rarray[i] = 3. - (10./(3.*L))*float(i)
if ((9.*L/10.) <= i < L):
rarray[i] = 0.
h5dset = h5grp.create_dataset('Data', data=rarray)
The Python Script¶
def thermal_transient_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.8, pg.21
# "A Transient 1D Temperature Distribution, Non-Zero Initial Temperature,
# Boundary Conditions of 1st and 2nd Kind"
# With some modification of the parameter values given in Kolditz (2015)
#
# Author: Jenn Frederick
# Date: 07/13/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
nxyz[0] = 10 # [m] nx
dxyz = lxyz/nxyz # [m]
K = 0.5787037 # [W/m-C]
Cp = 0.01 # [J/kg-C]
rho = 2000. # [kg/m^3]
chi = K/(Cp*rho)
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]
t_soln = np.array([0.015,0.05,0.10,0.50]) # [day]
T_soln = np.zeros((4,nx)) # [C]
T_pflotran = np.zeros((4,nx)) # [C]
# create the analytical solution
for time in range(4):
t = t_soln[time]*24.0*3600.0 # [sec]
sum_term_old = np.zeros(nx)
sum_term = np.zeros(nx)
n = 1
epsilon = 1.0
while epsilon > epsilon_value:
sum_term_old = sum_term
sum_term = sum_term_old + (np.cos(n*math.pi*x_soln/Lx)*np.exp(-chi*pow(n,2)*pow(math.pi,2)*t/pow(Lx,2))*(80./(3.*pow((n*math.pi),2)))*np.cos(n*math.pi/2.)*np.sin(n*math.pi/4.)*np.sin(3.*n*math.pi/20.))
epsilon = np.max(np.abs(sum_term_old-sum_term))
n = n + 1
T_soln[time,:] = 0.50 + sum_term
# run PFLOTRAN simulation
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)
index_string = 'Time: 1.50000E-02 d/Temperature [C]'
T_pflotran[0,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_1st_2nd_kind.h5',index_string,False)
index_string = 'Time: 5.00000E-02 d/Temperature [C]'
T_pflotran[1,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_1st_2nd_kind.h5',index_string,False)
index_string = 'Time: 1.00000E-01 d/Temperature [C]'
T_pflotran[2,:] = read_pflotran_output_1D(path+
'/1D_transient_thermal_BC_1st_2nd_kind.h5',index_string,False)
index_string = 'Time: 5.00000E-01 d/Temperature [C]'
T_pflotran[3,:] = read_pflotran_output_1D(path+
'/1D_transient_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.
dxyz = lxyz/nxyz
# Plot the PFLOTRAN and analytical solutions
plot_1D_transient(path,t_soln,'day',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 Transient Thermal Conduction, BCs of 1st & 2nd Kind¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
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.9, pg.23, “A Transient 2D Temperature Distribution, Non-Zero Initial Temperature, Boundary Conditions of 1st and 2nd Kind.”
The domain is a 100x100x1 meter square plate extending along the positive x and y axis and is made up of 50x50x1 hexahedral grid cells with dimensions 2x2x1 meters. The domain is composed of a single material and is assigned the following properties: thermal conductivity K = 0.5787037 W/(m-C); specific heat capacity Cp = 0.01 J/(m-C); density \(\rho\) = 2,000 kg/m^3.
The temperature is initially distributed according to
\[ \begin{align}\begin{aligned}T(t=0) = T_b f(x) f(y) + T_{offset}\\f(x) = \sum_{n=1}^{\infty} \left({80 \over {3(n\pi)^2}}\right)sin{{n\pi} \over 2} sin{{n \pi x} \over L} sin{{n\pi} \over 4} sin{{3n\pi} \over 20}\\f(y) = {1 \over 2} + \sum_{n=1}^{\infty} \left({80 \over {3(n\pi)^2}}\right) cos{{n \pi y} \over L} cos{{n\pi} \over 2} sin{{n\pi} \over 4} sin{{3n\pi} \over 20}\\\chi = {K \over {\rho c_p}}\end{aligned}\end{align} \]
where \(T_b\) = 1C and \(T_{offset}\) = 0.1 C. The boundary conditions are:
\[ \begin{align}\begin{aligned}T(0,y) = 0 \hspace{0.25in} x=0 \hspace{0.15in} face\\T(L,y) = 0 \hspace{0.25in} x=L \hspace{0.15in} face\\q(x,0) = 0 \hspace{0.25in} y=0 \hspace{0.15in} face\\q(x,L) = 0 \hspace{0.25in} y=L \hspace{0.15in} face\end{aligned}\end{align} \]
The transient temperature distribution is governed by,
\[\rho c_p {{\partial T} \over {\partial t}} = K \left({ {{\partial^{2} T} \over {\partial x^{2}}} + {{\partial^{2} T} \over {\partial y^{2}}} }\right)\]
The solution is given by,
\[ \begin{align}\begin{aligned}T(x,y,t) = T_b f(x,t) f(y,t) + T_{offset}\\f(x) = \sum_{n=1}^{\infty} exp\left({-\chi n^2 \pi^2 {t \over L^2}}\right)\left({80 \over {3(n\pi)^2}}\right) sin{{n \pi x} \over L} sin{{n\pi} \over 2} sin{{n\pi} \over 4} sin{{3n\pi} \over 20}\\f(y) = {1 \over 2} + \sum_{n=1}^{\infty} exp\left({-\chi n^2 \pi^2 {t \over L^2}}\right)\left({80 \over {3(n\pi)^2}}\right) cos{{n \pi y} \over L} cos{{n\pi} \over 2} sin{{n\pi} \over 4} sin{{3n\pi} \over 20}\\\chi = {K \over {\rho c_p}}\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 40 40 1
DXYZ
2.5d0
2.5d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 100.0d0 1.0d0
/
END
REGION west_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 100.0d0 1.0d0
/
END
REGION east_face
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 100.0d0 1.0d0
/
END
REGION north_face
FACE NORTH
COORDINATES
0.d0 100.0d0 0.d0
100.0d0 100.0d0 1.0d0
/
END
REGION south_face
FACE SOUTH
COORDINATES
0.d0 0.d0 0.d0
100.0d0 0.d0 1.0d0
/
END
MATERIAL_PROPERTY plate
ID 1
POROSITY 1.d-10
TORTUOSITY 1.d0
ROCK_DENSITY 2.0E3
HEAT_CAPACITY 0.01 J/kg-C
THERMAL_CONDUCTIVITY_DRY 0.5787037 W/m-C
THERMAL_CONDUCTIVITY_WET 0.5787037 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
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
#ENTHALPY CONSTANT 1.8890d0 MJ/kg
VISCOSITY CONSTANT 1.d-3 Pa-s
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
STRATA
REGION all
MATERIAL plate
END
TIME
FINAL_TIME 0.10 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.0001 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
TIMES day 0.04 0.06
#NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET temperature_initial
HDF5_DATASET_NAME initial
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE DATASET temperature_initial
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-20
END
FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 0.1d0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-20
END
FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 0.1d0 C
LIQUID_PRESSURE 101325 Pa
MOLE_FRACTION 1.d-20
END
FLOW_CONDITION north_face
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX 0.d0 W/m^2
LIQUID_FLUX 0.d0 m/yr
GAS_FLUX 0.d0 m/yr
END
FLOW_CONDITION south_face
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX 0.d0 W/m^2
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 40 40 1
DXYZ
2.5d0
2.5d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 100.0d0 1.0d0
/
END
REGION west_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 100.0d0 1.0d0
/
END
REGION east_face
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 100.0d0 1.0d0
/
END
REGION north_face
FACE NORTH
COORDINATES
0.d0 100.0d0 0.d0
100.0d0 100.0d0 1.0d0
/
END
REGION south_face
FACE SOUTH
COORDINATES
0.d0 0.d0 0.d0
100.0d0 0.d0 1.0d0
/
END
MATERIAL_PROPERTY plate
ID 1
POROSITY 1.d-10
TORTUOSITY 1.d0
ROCK_DENSITY 2.0E3
SPECIFIC_HEAT 0.01 J/kg-C
THERMAL_CONDUCTIVITY_DRY 0.5787037 W/m-C
THERMAL_CONDUCTIVITY_WET 0.5787037 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
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
#ENTHALPY CONSTANT 1.8890d0 MJ/kg
VISCOSITY CONSTANT 1.d-3 Pa-s
END
STRATA
REGION all
MATERIAL plate
END
TIME
FINAL_TIME 0.10 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.0001 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
TIMES day 0.04 0.06
#NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET temperature_initial
HDF5_DATASET_NAME initial
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE DATASET temperature_initial
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 0.1d0 C
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 0.1d0 C
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION north_face
TYPE
ENERGY_FLUX neumann
LIQUID_PRESSURE dirichlet
/
ENERGY_FLUX 0.d0 W/m^2
LIQUID_PRESSURE 101325 Pa
END
FLOW_CONDITION south_face
TYPE
ENERGY_FLUX neumann
LIQUID_PRESSURE dirichlet
/
ENERGY_FLUX 0.d0 W/m^2
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
ONLY_ENERGY_EQ
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
import math
import matplotlib.pyplot as plt
filename = 'dataset.h5'
h5file = File(filename,mode='w')
T_offset = 0.10 # [C]
T0 = 1. # [C]
L = 100. # [m]
dx = 1.0 # [m]
dy = 1.0 # [m]
dz = 1.0 # [m]
# 2D Surface:
# -----------------------------------------------------------------------------
# Temperature initial condition z=0 face
h5grp = h5file.create_group('initial')
h5grp.attrs['Dimension'] = np.bytes_('XY')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dx,dy]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
# Load the dataset values
nx = L*dx + 1
ny = L*dy + 1
rarray = np.zeros((nx,ny),'=f8')
fx = np.zeros(nx,'=f8')
for i in range(int(nx)):
if (0. <= i < (L/10.)):
fx[i] = 0.
if ((L/10.) <= i < (4.*L/10.)):
fx[i] = (10./(3.*L))*float(i) - (1./3.)
if ((4.*L/10.) <= i < (6.*L/10.)):
fx[i] = 1.
if ((6.*L/10.) <= i < (9.*L/10.)):
fx[i] = 3. - (10./(3.*L))*float(i)
if ((9.*L/10.) <= i < L):
fx[i] = 0.
fy = np.zeros(ny,'=f8')
for j in range(int(ny)):
if (0. <= j < (L/10.)):
fy[j] = 0.
if ((L/10.) <= j < (4.*L/10.)):
fy[j] = (10./(3.*L))*float(j) - (1./3.)
if ((4.*L/10.) <= j < (6.*L/10.)):
fy[j] = 1.
if ((6.*L/10.) <= j < (9.*L/10.)):
fy[j] = 3. - (10./(3.*L))*float(j)
if ((9.*L/10.) <= j < L):
fy[j] = 0.
for i in range(int(nx)):
for j in range(int(ny)):
rarray[i][j] = T0*fx[i]*fy[j] + T_offset
h5dset = h5grp.create_dataset('Data', data=rarray)
The Python Script¶
def thermal_transient_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.9, pg.23
# "A Transient 2D Temperature Distribution, Non-Zero Initial Temperature,
# Boundary Conditions of 1st and 2nd Kind"
#
# Author: Jenn Frederick
# Date: 07/13/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] = 100. # [m] ly
nxyz[0] = 10 # [m] nx
nxyz[1] = 10 # [m] ny
dxyz = lxyz/nxyz # [m]
T0 = 1. # [C]
T_offset = 0.1 # [C]
K = 0.5787037 # [W/m-C]
rho = 2000. # [kg/m^3]
Cp = 0.01 # [J/kg-C]
chi = K/(rho*Cp)
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((4,ny,nx)) # [C]
T_pflotran = np.zeros((4,ny,nx)) # [C]
t_soln = np.array([0.00,0.04,0.06,0.10]) # [day]
x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx) # [m]
y_soln = np.linspace(0.+(dy/2.),Ly-(dy/2.),ny) # [m]
T1x = np.zeros(int(nx))
T2y = np.zeros(int(ny))
# for time=0, use the actual initial condition as the analytical
# solution, not the analytical solution with t=0:
fx = np.zeros(int(nx))
for i in range(int(nx)):
x = x_soln[i]
if (0. <= x < (Lx/10.)):
fx[i] = 0.
if ((Lx/10.) <= x < (4.*Lx/10.)):
fx[i] = (10./(3.*Lx))*(x) - (1./3.)
if ((4.*Lx/10.) <= x < (6.*Lx/10.)):
fx[i] = 1.
if ((6.*Lx/10.) <= x < (9.*Lx/10.)):
fx[i] = 3. - (10./(3.*Lx))*(x)
if ((9.*Lx/10.) <= x < Lx):
fx[i] = 0.
fy = np.zeros(int(ny))
for j in range(int(ny)):
y = y_soln[j]
if (0. <= y < (Ly/10.)):
fy[j] = 0.
if ((Ly/10.) <= y < (4.*Ly/10.)):
fy[j] = (10./(3.*Ly))*(y) - (1./3.)
if ((4.*Ly/10.) <= y < (6.*Ly/10.)):
fy[j] = 1.
if ((6.*Ly/10.) <= y < (9.*Ly/10.)):
fy[j] = 3. - (10./(3.*Ly))*(y)
if ((9.*Ly/10.) <= y < Ly):
fy[j] = 0.
T_soln_t0 = np.zeros((nx,ny))
for i in range(int(nx)):
for j in range(int(ny)):
T_soln_t0[i][j] = T0*fx[i]*fy[j] + T_offset
# create the analytical solution
for time in range(4):
t = t_soln[time]*24.0*3600.0 # [sec]
# create T1y
sum_term_y = np.zeros(int(ny))
sum_term_old_y = np.zeros(int(ny))
n = 1
epsilon = 1
while epsilon > epsilon_value:
sum_term_old_y = sum_term_y
sum_term_y = sum_term_old_y + (np.cos(n*math.pi*y_soln/Ly)*np.exp(-chi*pow(n,2)*pow(math.pi,2)*t/pow(Ly,2))*(80./(3.*pow((n*math.pi),2)))*np.cos(n*math.pi/2.)*np.sin(n*math.pi/4.)*np.sin(3.*n*math.pi/20.))
epsilon = np.max(np.abs(sum_term_old_y-sum_term_y))
n = n + 1
T2y = 0.5 + sum_term_y
# create T1x
sum_term_x = np.zeros(int(nx))
sum_term_old_x = np.zeros(int(nx))
n = 1
epsilon = 1
while epsilon > epsilon_value:
sum_term_old_x = sum_term_x
sum_term_x = sum_term_old_x + (np.sin(n*math.pi*x_soln/Lx)*np.exp(-chi*pow(n,2)*pow(math.pi,2)*t/pow(Lx,2))*(80./(3.*pow((n*math.pi),2)))*np.sin(n*math.pi/2.)*np.sin(n*math.pi/4.)*np.sin(3.*n*math.pi/20.))
epsilon = np.max(np.abs(sum_term_old_x-sum_term_x))
n = n + 1
T1x = sum_term_x
for i in range(int(nx)):
for j in range(int(ny)):
T_soln[time,i,j] = T0*T1x[i]*T2y[j] + T_offset
T_soln[0,:,:] = T_soln_t0[:,:]
# run PFLOTRAN simulation
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)
index_string = 'Time: 0.00000E+00 d/Temperature [C]'
T_pflotran[0,:,:] = read_pflotran_output_2D(path+
'/2D_transient_thermal_BC_1st_2nd_kind.h5',index_string,False)
index_string = 'Time: 4.00000E-02 d/Temperature [C]'
T_pflotran[1,:,:] = read_pflotran_output_2D(path+
'/2D_transient_thermal_BC_1st_2nd_kind.h5',index_string,False)
index_string = 'Time: 6.00000E-02 d/Temperature [C]'
T_pflotran[2,:,:] = read_pflotran_output_2D(path+
'/2D_transient_thermal_BC_1st_2nd_kind.h5',index_string,False)
index_string = 'Time: 1.00000E-01 d/Temperature [C]'
T_pflotran[3,:,:] = read_pflotran_output_2D(path+
'/2D_transient_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 solution:
plot_2D_transient(path,t_soln,'day',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.
