1D Steady Flow (Pressure), BCs of 1st Kind¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
The PFLOTRAN Input File (RICHARDS 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.2.1, pg.27, “A 1D Steady-State Pressure 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 materials are assigned the following properties: porosity = 0.5; fluid viscosity \(\mu\) = 1.0 mPa-s; permeability = 1e-15 m^2; fluid density rho = 1,000 kg/m^3.
The pressure is initially uniform at p(t=0) = 1.5 MPa. The boundary pressures are p(x=0) = 2.0 MPa and p(x=L) = 1.0 MPa, where L = 100 m. The simulation is run until the steady-state pressure distribution develops.
The LaPlace equation governs the steady-state pressure distribution,
\[{{\partial^{2} p} \over {\partial x^{2}}} = 0\]
The solution is given by,
\[p(x) = ax + b\]
When the boundary conditions p(x=0) = 2.0 MPa and p(x=L) = 1.0 MPa are applied, the solution becomes,
\[p(x) = {-x \over L} + 2.0\]
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.¶
Comparison of the PFLOTRAN vs. analytical solution for RICHARDS 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
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 0.50d0
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-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.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
STRATA
REGION all
MATERIAL beam
END
TIME
FINAL_TIME 10 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 1.d0 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.5 MPa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION left_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 2 MPa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION right_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1 MPa
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 10 1 1
DXYZ
10.0d0
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 0.50d0
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-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.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
STRATA
REGION all
MATERIAL beam
END
TIME
FINAL_TIME 10 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 1.d0 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.5 MPa
END
FLOW_CONDITION left_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 2 MPa
END
FLOW_CONDITION right_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1 MPa
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 (RICHARDS Mode)¶
The RICHARDS Mode PFLOTRAN input file can be downloaded
here.
SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE RICHARDS
/
/
END
SUBSURFACE
GRID
TYPE structured
NXYZ 10 1 1
DXYZ
10.0d0
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 0.50d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
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.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
STRATA
REGION all
MATERIAL beam
END
TIME
FINAL_TIME 10 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 1.d0 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 1.5 MPa
END
FLOW_CONDITION left_end
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 2 MPa
END
FLOW_CONDITION right_end
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 1 MPa
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 flow_steady_1D_BC1stkind(path,input_prefix,remove,screen_on,pf_exe,
mpi_exe,num_tries):
#==============================================================================#
#
# 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.2.1, pg.27
# "A 1D Steady-State Pressure Distribution, Boundary Conditions of 1st Kind"
#
# Author: Jenn Frederick
# Date: 07/21/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]
P0 = 2. # [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]
p_soln = np.zeros(nx+2) # [MPa]
p_pflotran = np.zeros(nx) # [Pa]
# create the analytical solution
p_soln = np.array(-x_soln/Lx + P0) # [MPa]
# 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 d/Liquid_Pressure [Pa]'
p_pflotran[:] = read_pflotran_output_1D(path+
'/1D_steady_pressure_BC_1st_kind.h5',index_string,remove)
ierr = check(p_pflotran)
# convert pressure units [Pa -> MPa]
p_pflotran = p_pflotran/1.0e6 # [MPa]
max_percent_error = calc_relative_error(p_soln[1:nx+1],p_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,p_soln,x_pflotran,p_pflotran,'Distance [m]',
'Pressure [MPa]',"{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 Flow (Pressure), BCs of 1st & 2nd Kind¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
The PFLOTRAN Input File (RICHARDS 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.2.2, pg.27, “A 1D Steady-State Pressure 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 20x1x1 cubic grid cells with dimensions 5x10x10 meters. The domain is composed of two materials and is assigned the following properties: permeability k1(x<2L/5) = 1e-12 m^2; permeability k2(x>2L/5) = 3e-12 m^2; specific heat capacity Cp = 0.001 J/(m-C); rock density = 2,800 kg/m^3; fluid density = 1,000 kg/m^3; fluid viscosity \(\mu\) = 1 mPa-s; porosity = 0.50.
The temperature is initially uniform at T(t=0) = 25C, and the pressure is initially uniform at p(t=0) = 1.0 MPa. The temperature at both ends is T(x=0,x=L) = 25.0C. The pressure at the left end is p(x=0) = 1.0 MPa, and a fluid flux of q(x=L) = -1.5e-5 m/s is applied at the right end, where L = 100 m. The simulation is run until the steady-state pressure distribution develops.
The LaPlace equation governs the steady-state pressure distribution,
\[{{\partial^{2} p} \over {\partial x^{2}}} = 0\]
The solution is given by,
\[ \begin{align}\begin{aligned}p(x<2L/5) = a_1 x + b_1\\p(x>2L/5) = a_2 x + b_2\end{aligned}\end{align} \]
When the boundary conditions p(x=0) = 1.0 MPa and q(x=L) = -1.5e-5 m/s are applied, the solution becomes,
\[ \begin{align}\begin{aligned}p(x<2L/5) = -{{q \mu x} \over {k1}} + p(x=0)\\p(x>2L/5) = -{{q \mu x} \over {k2}} + p(x=0) + {{2 q \mu L} \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.¶
Comparison of the PFLOTRAN vs. analytical solution for RICHARDS 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
1.0d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 1.0d0 1.0d0
/
END
REGION left
COORDINATES
0.d0 0.d0 0.d0
40.0d0 1.0d0 1.0d0
/
END
REGION right
COORDINATES
40.0d0 0.0d0 0.0d0
100.0d0 1.0d0 1.0d0
/
END
REGION left_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END
REGION right_face
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 1.0d0 1.0d0
/
END
MATERIAL_PROPERTY beam_left
ID 1
POROSITY 0.5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3 kg/m^3
HEAT_CAPACITY 1.d-3 J/kg-C
THERMAL_CONDUCTIVITY_DRY 100 W/m-C
THERMAL_CONDUCTIVITY_WET 100 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-12
PERM_Y 1.d-12
PERM_Z 1.d-12
/
END
MATERIAL_PROPERTY beam_right
ID 2
POROSITY 0.5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3 kg/m^3
HEAT_CAPACITY 1.d-3 J/kg-C
THERMAL_CONDUCTIVITY_DRY 100 W/m-C
THERMAL_CONDUCTIVITY_WET 100 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 3.d-12
PERM_Y 3.d-12
PERM_Z 3.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 left
MATERIAL beam_left
END
STRATA
REGION right
MATERIAL beam_right
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 10 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.25 day at 0.d0 day
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 25.d0 C
LIQUID_PRESSURE 1.0 MPa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION left_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.0 MPa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION right_face
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX 0.d0 W/m^2
LIQUID_FLUX +1.50d-5 m/sec
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
1.0d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 1.0d0 1.0d0
/
END
REGION left
COORDINATES
0.d0 0.d0 0.d0
40.0d0 1.0d0 1.0d0
/
END
REGION right
COORDINATES
40.0d0 0.0d0 0.0d0
100.0d0 1.0d0 1.0d0
/
END
REGION left_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END
REGION right_face
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 1.0d0 1.0d0
/
END
MATERIAL_PROPERTY beam_left
ID 1
POROSITY 0.5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3 kg/m^3
SPECIFIC_HEAT 1.d-3 J/kg-C
THERMAL_CONDUCTIVITY_DRY 100 W/m-C
THERMAL_CONDUCTIVITY_WET 100 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-12
PERM_Y 1.d-12
PERM_Z 1.d-12
/
END
MATERIAL_PROPERTY beam_right
ID 2
POROSITY 0.5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3 kg/m^3
SPECIFIC_HEAT 1.d-3 J/kg-C
THERMAL_CONDUCTIVITY_DRY 100 W/m-C
THERMAL_CONDUCTIVITY_WET 100 W/m-C
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 3.d-12
PERM_Y 3.d-12
PERM_Z 3.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 left
MATERIAL beam_left
END
STRATA
REGION right
MATERIAL beam_right
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 10 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.25 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.0 MPa
END
FLOW_CONDITION left_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.0 MPa
END
FLOW_CONDITION right_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX +1.50d-5 m/sec
#LIQUID_FLUX +1.7d-5 m/sec matches!
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 (RICHARDS Mode)¶
The RICHARDS Mode PFLOTRAN input file can be downloaded
here.
SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE RICHARDS
/
/
END
SUBSURFACE
GRID
TYPE structured
NXYZ 10 1 1
DXYZ
10.0d0
1.0d0
1.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
100.0d0 1.0d0 1.0d0
/
END
REGION left
COORDINATES
0.d0 0.d0 0.d0
40.0d0 1.0d0 1.0d0
/
END
REGION right
COORDINATES
40.0d0 0.0d0 0.0d0
100.0d0 1.0d0 1.0d0
/
END
REGION left_face
FACE WEST
COORDINATES
0.d0 0.d0 0.d0
0.d0 1.0d0 1.0d0
/
END
REGION right_face
FACE EAST
COORDINATES
100.0d0 0.d0 0.d0
100.0d0 1.0d0 1.0d0
/
END
MATERIAL_PROPERTY beam_left
ID 1
POROSITY 0.5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3 kg/m^3
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-12
PERM_Y 1.d-12
PERM_Z 1.d-12
/
END
MATERIAL_PROPERTY beam_right
ID 2
POROSITY 0.5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3 kg/m^3
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 3.d-12
PERM_Y 3.d-12
PERM_Z 3.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 left
MATERIAL beam_left
END
STRATA
REGION right
MATERIAL beam_right
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 10 day
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.25 day at 0.d0 day
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 1.0 MPa
END
FLOW_CONDITION left_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 1.0 MPa
END
FLOW_CONDITION right_face
TYPE
LIQUID_FLUX neumann
/
LIQUID_FLUX +1.50d-5 m/sec
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 flow_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.2.2, pg.27
# "A 1D Steady-State Pressure Distribution, Boundary Conditions of 1st and
# 2nd Kind"
#
# Author: Jenn Frederick
# Date: 07/21/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]
K1 = 1.e-12 # [m^2]
K2 = 3.e-12 # [m^2]
mu = 1.e-3 # [Pa-s]
q = -1.5e-5 # [m/s]
p_at0 = 1.e6 # [Pa]
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]
p_soln = np.zeros(nx+1) # [Pa]
p_pflotran = np.zeros(nx) # [Pa]
# create the analytical solution
k = -1
for x in x_soln:
k = k + 1
if x <= (2.*Lx/5.):
p_soln[k] = np.array( -(((q*mu)/K1)*x) + p_at0 )
else:
p_soln[k] = np.array( -(((q*mu)/K2)*x) + p_at0 + (q*mu*(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+01 d/Liquid_Pressure [Pa]'
p_pflotran[:] = read_pflotran_output_1D(path+
'/1D_steady_pressure_BC_1st_2nd_kind.h5',index_string,remove)
ierr = check(p_pflotran)
# Convert pressure units [Pa -> MPa]
p_pflotran = p_pflotran/1.e6 # [MPa]
p_soln = p_soln/1.e6 # [MPa]
max_percent_error = calc_relative_error(p_soln[1:nx+1],p_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,p_soln,x_pflotran,p_pflotran,'Distance [m]',
'Pressure [MPa]',"{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 Flow (Pressure), Hydrostatic¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
The PFLOTRAN Input File (RICHARDS 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.2.6, pg.32, “A Hydrostatic Pressure Distribution.”
The domain is a 100x20x20 meter rectangular beam extending along the positive z-axis and is made up of 1x1x100 hexahedral grid cells with dimensions 20x20x1 meters. The material properties should not influence the solution for this test problem, except for the fluid density. The following properties are assigned: porosity = 0.5; fluid viscosity \(\mu\) = 1.0 mPa-s; permeability = 1e-12 m^2; fluid density rho = 1,019 kg/m^3.
The pressure is initially uniform at p(t=0) = 0.5 MPa. The top pressure boundary condition is p(z=H) = 0 MPa and the bottom of the column is assigned a no flow boundary condition, q(z=0) = 0 m/s, where H = 100 m. The simulation is run until the steady-state pressure distribution develops.
The solution is given by the hydrostatic pressure distribution, and does not depend on the material properties assigned,
\[p(z) = \rho g (H - z)\]
where g is the standard value for gravity (9.81 m/s2).
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.¶
Comparison of the PFLOTRAN vs. analytical solution for RICHARDS 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 1 1 10
DXYZ
1.0d0
1.0d0
10.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
1.0d0 1.0d0 100.0d0
/
END
REGION top_end
FACE TOP
COORDINATES
0.d0 0.d0 100.0d0
1.0d0 1.0d0 100.0d0
/
END
REGION bottom_end
FACE BOTTOM
COORDINATES
0.d0 0.d0 0.d0
1.0d0 1.0d0 0.d0
/
END
MATERIAL_PROPERTY beam
ID 1
POROSITY 0.50d0
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 beam
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1019.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.1d0 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 25.d0 C
LIQUID_PRESSURE 101325. Pa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION top_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 101325. Pa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION bottom_end
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0. m/sec
GAS_FLUX 0. m/sec
END
INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END
BOUNDARY_CONDITION top_end
REGION top_end
FLOW_CONDITION top_end
END
BOUNDARY_CONDITION bottom_end
FLOW_CONDITION bottom_end
REGION bottom_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 1 1 10
DXYZ
1.0d0
1.0d0
10.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
1.0d0 1.0d0 100.0d0
/
END
REGION top_end
FACE TOP
COORDINATES
0.d0 0.d0 100.0d0
1.0d0 1.0d0 100.0d0
/
END
REGION bottom_end
FACE BOTTOM
COORDINATES
0.d0 0.d0 0.d0
1.0d0 1.0d0 0.d0
/
END
MATERIAL_PROPERTY beam
ID 1
POROSITY 0.50d0
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 beam
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1019.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.1d0 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 25.d0 C
LIQUID_PRESSURE 101325. Pa
END
FLOW_CONDITION top_end
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 101325. Pa
END
FLOW_CONDITION bottom_end
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0 m/sec
END
INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END
BOUNDARY_CONDITION top_end
REGION top_end
FLOW_CONDITION top_end
END
BOUNDARY_CONDITION bottom_end
FLOW_CONDITION bottom_end
REGION bottom_end
END
END_SUBSURFACE
The PFLOTRAN Input File (RICHARDS Mode)¶
The RICHARDS Mode PFLOTRAN input file can be downloaded
here.
SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE RICHARDS
/
/
END
SUBSURFACE
GRID
TYPE structured
NXYZ 1 1 10
DXYZ
1.0d0
1.0d0
10.0d0
END
END
REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
1.0d0 1.0d0 100.0d0
/
END
REGION top_end
FACE TOP
COORDINATES
0.d0 0.d0 100.0d0
1.0d0 1.0d0 100.0d0
/
END
REGION bottom_end
FACE BOTTOM
COORDINATES
0.d0 0.d0 0.d0
1.0d0 1.0d0 0.d0
/
END
MATERIAL_PROPERTY beam
ID 1
POROSITY 0.50d0
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
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 beam
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1019.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.1d0 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END
FLOW_CONDITION initial
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 101325. Pa
END
FLOW_CONDITION top_end
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 101325. Pa
END
FLOW_CONDITION bottom_end
TYPE
LIQUID_FLUX neumann
/
LIQUID_FLUX 0. m/sec
END
INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END
BOUNDARY_CONDITION top_end
REGION top_end
FLOW_CONDITION top_end
END
BOUNDARY_CONDITION bottom_end
FLOW_CONDITION bottom_end
REGION bottom_end
END
END_SUBSURFACE
The Python Script¶
def flow_steady_1D_hydrostatic(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.2.6, pg.32
# "A Hydrostatic Pressure Distribution"
#
# Author: Jenn Frederick
# Date: 08/23/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[2] = 100. # [m] lz
nxyz[2] = 10 # [m] nz
dxyz = lxyz/nxyz # [m]
g = 9.81 # [m/s^2]
rho = 1019.0 # [kg/m^3]
p_offset = 101325.0 # [Pa]
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
z_soln = np.linspace(0.+(dz/2.),Lz-(dz/2.),nz) # [m]
z_pflotran = z_soln # [m]
z_soln = np.concatenate(([0.],z_soln),axis=0) # [m]
p_soln = np.zeros(nz+1) # [Pa]
p_pflotran = np.zeros((nx,ny,nz)) # [Pa] needs to be 3D
# create the analytical solution
k = -1
for z in z_soln:
k = k + 1
p_soln[k] = rho*g*(Lz-z_soln[k]) + p_offset
# 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+00 y/Liquid_Pressure [Pa]'
p_pflotran[:,:,:] = read_pflotran_output_3D(path+
'/1D_steady_pressure_hydrostatic.h5',index_string,remove)
ierr = check(p_pflotran)
# Convert pressure units [Pa -> MPa]
p_pflotran = p_pflotran/1.e6 # [MPa]
p_soln = p_soln/1.e6 # [MPa]
max_percent_error = calc_relative_error(p_soln[1:nz+1],p_pflotran[0,0,:],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[2] = nxyz[2]*2.
dxyz = lxyz/nxyz
# Plot the PFLOTRAN and analytical solutions
plot_1D_steady(path,z_soln,p_soln,z_pflotran,p_pflotran[0,0,:],'Distance [m]',
'Pressure [MPa]',"{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 Flow (Pressure), BCs of 1st Kind¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
The PFLOTRAN Input File (RICHARDS 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.2.3, pg.29, “A 2D Steady-State Pressure 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 20x20x1 hexahedral grid cells with dimensions 0.05x0.05x1 meters. The materials are assigned the following properties: porosity = 0.5; permeability = 1e-15 m2; fluid density rho = 1000 kg/m3; fluid viscosity \(\mu\) = 1 mPa-s.
The pressure is initially uniform at p(t=0) = 1.0 MPa. The boundary pressure conditions (in units of MPa) are:
\[ \begin{align}\begin{aligned}p(0,y) = 1 \hspace{0.25in} x=0 \hspace{0.15in} face\\p(L,y) = p_0 {y \over L} + 1 \hspace{0.25in} x=L \hspace{0.15in} face\\p(x,0) = 1 \hspace{0.25in} y=0 \hspace{0.15in} face\\p(x,L) = p_0 {x \over L} + 1 \hspace{0.25in} y=L \hspace{0.15in} face\end{aligned}\end{align} \]
where L = 1 m and p0 = 2.0 MPa. The simulation is run until the steady-state pressure distribution develops.
The LaPlace equation governs the steady-state pressure distribution,
\[{{\partial^{2} p} \over {\partial x^{2}}} + {{\partial^{2} p} \over {\partial y^{2}}}= 0\]
After applying these boundary conditions, the solution is given by,
\[p(x,y) = {x \over L}{y \over L} + p_0\]
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.¶
Comparison of the PFLOTRAN vs. analytical solution for RICHARDS 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.5
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-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
STRATA
REGION all
MATERIAL plate
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.1 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END
DATASET pressure_bc_x
HDF5_DATASET_NAME x_line_node_centered
FILENAME ../dataset.h5
END
DATASET pressure_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 25.d0 C
LIQUID_PRESSURE 1.d0 MPa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.d0 MPa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_y
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_x
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.d0 MPa
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 south_face
REGION south_face
FLOW_CONDITION south_face
END
BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END
BOUNDARY_CONDITION north_face
FLOW_CONDITION north_face
REGION north_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.5
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-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
STRATA
REGION all
MATERIAL plate
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.1 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END
DATASET pressure_bc_x
HDF5_DATASET_NAME x_line_node_centered
FILENAME ../dataset.h5
END
DATASET pressure_bc_y
HDF5_DATASET_NAME y_line_node_centered
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.d0 MPa
END
FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.d0 MPa
END
FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_y
END
FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_x
END
FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.d0 MPa
END
INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END
BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END
BOUNDARY_CONDITION south_face
REGION south_face
FLOW_CONDITION south_face
END
BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END
BOUNDARY_CONDITION north_face
FLOW_CONDITION north_face
REGION north_face
END
END_SUBSURFACE
The PFLOTRAN Input File (RICHARDS Mode)¶
The Richards Mode PFLOTRAN input file can be downloaded
here.
SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE RICHARDS
/
/
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.5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
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
STRATA
REGION all
MATERIAL plate
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.1 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END
DATASET pressure_bc_x
HDF5_DATASET_NAME x_line_node_centered
FILENAME ../dataset.h5
END
DATASET pressure_bc_y
HDF5_DATASET_NAME y_line_node_centered
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 1.d0 MPa
END
FLOW_CONDITION west_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 1.d0 MPa
END
FLOW_CONDITION east_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_y
END
FLOW_CONDITION north_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_x
END
FLOW_CONDITION south_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 1.d0 MPa
END
INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END
BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END
BOUNDARY_CONDITION south_face
REGION south_face
FLOW_CONDITION south_face
END
BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END
BOUNDARY_CONDITION north_face
FLOW_CONDITION north_face
REGION north_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')
Po = 1.0e6 # [Pa]
p_offset = 1.0e6 # [Pa]
L = 1.0 # [m]
# 1d line in x
# Pressure 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):
x = (float(i)*L)/(nx-1)
rarray[i] = Po*(x/L) + p_offset
h5dset = h5grp.create_dataset('Data', data=rarray)
# 1d line in y
# Pressure 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 j in range(ny):
y = (float(j)*L)/(ny-1)
rarray[j] = Po*(y/L) + p_offset
h5dset = h5grp.create_dataset('Data', data=rarray)
The Python Script¶
def flow_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.2.3, pg.29
# "A 2D Steady-State Pressure Distribution, Boundary Conditions of 1st Kind"
#
# Author: Jenn Frederick
# Date: 07/22/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] = 1. # [m] lx
lxyz[1] = 1. # [m] ly
nxyz[0] = 10 # [m] nx
nxyz[1] = 10 # [m] ny
dxyz = lxyz/nxyz # [m]
p0 = 1. # [MPa]
p_offset = 1.0 # [MPa]
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
p_soln = np.zeros((nx,ny)) # [MPa]
p_pflotran = np.zeros((nx,ny)) # [Pa]
x_soln = np.linspace(0.+(dx/2.),Lx-(dx/2.),nx) # [m]
y_soln = np.linspace(0.+(dy/2.),Ly-(dy/2.),ny) # [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
p_soln[i,j] = p0*(x/Lx)*(y/Ly) + p_offset # [MPa]
# 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+00 y/Liquid_Pressure [Pa]'
p_pflotran[:,:] = read_pflotran_output_2D(path+
'/2D_steady_pressure_BC_1st_kind.h5',index_string,remove)
ierr = check(p_pflotran)
# Convert pressure units [Pa -> MPa]
p_pflotran = p_pflotran/1.e6 # [MPa]
max_percent_error = calc_relative_error(p_soln,p_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,p_soln,p_pflotran,'Distance [m]',
'Pressure [MPa]',"{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 Flow (Pressure), BCs of 1st & 2nd Kind¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
The PFLOTRAN Input File (RICHARDS 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.2.4, pg.29, “A 2D Steady-State Pressure Distribution, Boundary Conditions of 1st and 2nd Kind.”
The domain is a 2x1x1 meter rectangular slab extending along the positive x- 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: permeability k = 1e-14 m^2; rock density = 2,800 kg/m^3; fluid density = 1,000 kg/m^3; fluid viscosity \(\mu\) = 1 mPa-s; porosity = 0.50.
The pressure is initially uniform at p(t=0) = 1.0 MPa. The pressure boundary conditions (in units of MPa) and fluid flux boundary conditions (in units of m/s) are:
\[ \begin{align}\begin{aligned}p(x,0) = p_0 {x \over L} \hspace{0.25in} y=0 \hspace{0.15in} face\\p(x,L) = p_0 {{x+2L} \over L} \hspace{0.25in} y=L \hspace{0.15in} face\\p(2L,y) = p_0 {{2y+2L} \over L} \hspace{0.25in} x=2L \hspace{0.15in} face\\q(0,y) = {{-k}\over{\mu}}{p_0 \over L} = -10^{-5} \hspace{0.25in} x=0 \hspace{0.15in} face\end{aligned}\end{align} \]
where L = 1 m and p0 = 1 MPa. The simulation is run until the steady-state pressure distribution develops.
The LaPlace equation governs the steady-state pressure distribution,
\[{{\partial^{2} p} \over {\partial x^{2}}} + {{\partial^{2} p} \over {\partial y^{2}}} = 0\]
When the boundary conditions are applied, the solution becomes,
\[p(x,y) = {p_0 \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.¶
Comparison of the PFLOTRAN vs. analytical solution for RICHARDS 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 0.5
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-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 plate
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.1 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET pressure_bc_north
HDF5_DATASET_NAME x_line_node_centered_north
FILENAME ../dataset.h5
END
DATASET pressure_bc_south
HDF5_DATASET_NAME x_line_node_centered_south
FILENAME ../dataset.h5
END
DATASET pressure_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 10.d0 C
LIQUID_PRESSURE 3.d0 MPa
MOLE_FRACTION 1.d-20
END
FLOW_CONDITION west_face
TYPE
ENERGY_FLUX neumann
LIQUID_FLUX neumann
GAS_FLUX neumann
/
ENERGY_FLUX 0.d0 W/m^2
LIQUID_FLUX -1.d-5 m/s
GAS_FLUX 0.d0 m/yr
END
FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 10.d0 C
LIQUID_PRESSURE DATASET pressure_bc_east
MOLE_FRACTION 1.d-20
END
FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 10.d0 C
LIQUID_PRESSURE DATASET pressure_bc_north
MOLE_FRACTION 1.d-20
END
FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 10.d0 C
LIQUID_PRESSURE DATASET pressure_bc_south
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 south_face
REGION south_face
FLOW_CONDITION south_face
END
BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END
BOUNDARY_CONDITION north_face
FLOW_CONDITION north_face
REGION north_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.5
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 plate
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.1 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET pressure_bc_north
HDF5_DATASET_NAME x_line_node_centered_north
FILENAME ../dataset.h5
END
DATASET pressure_bc_south
HDF5_DATASET_NAME x_line_node_centered_south
FILENAME ../dataset.h5
END
DATASET pressure_bc_east
HDF5_DATASET_NAME y_line_node_centered_east
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 10.d0 C
LIQUID_PRESSURE 3.d0 MPa
END
FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
/
TEMPERATURE 10.d0 C
LIQUID_FLUX -1.d-5 m/s
END
FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 10.d0 C
LIQUID_PRESSURE DATASET pressure_bc_east
END
FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 10.d0 C
LIQUID_PRESSURE DATASET pressure_bc_north
END
FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 10.d0 C
LIQUID_PRESSURE DATASET pressure_bc_south
END
INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END
BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END
BOUNDARY_CONDITION south_face
REGION south_face
FLOW_CONDITION south_face
END
BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END
BOUNDARY_CONDITION north_face
FLOW_CONDITION north_face
REGION north_face
END
END_SUBSURFACE
The PFLOTRAN Input File (RICHARDS Mode)¶
The RICHARDS Mode PFLOTRAN input file can be downloaded
here.
SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE RICHARDS
/
/
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.5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
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 plate
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.1 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET pressure_bc_north
HDF5_DATASET_NAME x_line_node_centered_north
FILENAME ../dataset.h5
END
DATASET pressure_bc_south
HDF5_DATASET_NAME x_line_node_centered_south
FILENAME ../dataset.h5
END
DATASET pressure_bc_east
HDF5_DATASET_NAME y_line_node_centered_east
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 3.d0 MPa
END
FLOW_CONDITION west_face
TYPE
LIQUID_FLUX neumann
/
LIQUID_FLUX -1.d-5 m/s
END
FLOW_CONDITION east_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_east
END
FLOW_CONDITION north_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_north
END
FLOW_CONDITION south_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_south
END
INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END
BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION west_face
END
BOUNDARY_CONDITION south_face
REGION south_face
FLOW_CONDITION south_face
END
BOUNDARY_CONDITION east_face
FLOW_CONDITION east_face
REGION east_face
END
BOUNDARY_CONDITION north_face
FLOW_CONDITION north_face
REGION north_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.0 # [m]
p0 = 1.0e6 # [Pa]
p_offset = 1.0e6 # [Pa]
# 1d line in x
# Pressure 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):
x = (float(i)*2*L)/(nx-1)
rarray[i] = p0*x/L + p_offset
h5dset = h5grp.create_dataset('Data', data=rarray)
# 1d line in x
# Pressure 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):
x = (float(i)*2*L)/(nx-1)
rarray[i] = (p0/L)*(x+(2*L)) + p_offset
h5dset = h5grp.create_dataset('Data', data=rarray)
# 1d line in y
# Pressure 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):
y = (float(j)*L)/(ny-1)
rarray[j] = (p0/L)*((2*L)+(2*y)) + p_offset
h5dset = h5grp.create_dataset('Data', data=rarray)
The Python Script¶
def flow_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.2.4, pg.29
# "A 2D Steady-State Pressure Distribution, Boundary Conditions of 1st and
# 2nd Kind"
#
# Author: Jenn Frederick
# Date: 07/25/2017
# ******************************************************************************
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] = 2. # [m] lx
lxyz[1] = 1. # [m] ly
nxyz[0] = 20 # [m] nx
nxyz[1] = 10 # [m] ny
dxyz = lxyz/nxyz # [m]
p0 = 1. # [MPa]
p_offset = 1.0 # [MPa]
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
p_soln = np.zeros((nx,ny)) # [MPa]
p_pflotran = np.zeros((nx,ny)) # [Pa]
x_soln = np.linspace(0.+(dx/2.),(Lx)-(dx/2.),nx) # [m]
y_soln = np.linspace(0.+(dy/2.),Ly-(dy/2.),ny) # [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
p_soln[i,j] = (p0/Ly)*(x+(2*y)) + p_offset
# 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+00 y/Liquid_Pressure [Pa]'
p_pflotran[:,:] = read_pflotran_output_2D(path+
'/2D_steady_pressure_BC_1st_2nd_kind.h5',index_string,remove)
ierr = check(p_pflotran)
# Convert pressure units [Pa -> MPa]
p_pflotran = p_pflotran/1.e6 # [MPa]
max_percent_error = calc_relative_error(p_soln,p_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,p_soln,p_pflotran,'Distance [m]',
'Pressure [MPa]',"{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 Flow (Pressure), BCs of 1st Kind¶
The PFLOTRAN Input File (GENERAL Mode)
The PFLOTRAN Input File (TH Mode)
The PFLOTRAN Input File (RICHARDS 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.2.5, pg.31, “A 3D Steady-State Pressure 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: permeability k = 1e-10 m^2; rock density = 2,800 kg/m^3; fluid density = 1,000 kg/m^3; fluid viscosity \(\mu\) = 1 mPa-s; porosity = 0.50.
The pressure is initially uniform at p(t=0) = 1.0 MPa. The pressure boundary conditions (in units of MPa) are:
\[ \begin{align}\begin{aligned}p(0,y,z) = p_0 \left( {{0}+{y \over L}+{z \over L}} \right) \hspace{0.25in} x=0 \hspace{0.15in} face\\p(x,0,z) = p_0 \left( {{x \over L}+{0}+{z \over L}} \right) \hspace{0.25in} y=0 \hspace{0.15in} face\\
p(x,y,0) = p_0 \left( {{x \over L}+{y \over L}+{0}} \right) \hspace{0.25in} z=0 \hspace{0.15in} face\\
p(L,y,z) = p_0 \left( {{L}+{y \over L}+{z \over L}} \right) \hspace{0.25in} x=L \hspace{0.15in} face\\
p(x,L,z) = p_0 \left( {{x \over L}+{L}+{z \over L}} \right) \hspace{0.25in} y=L \hspace{0.15in} face\\
p(x,y,L) = p_0 \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 p0 = 1 MPa. The simulation is run until the steady-state pressure distribution develops.
The LaPlace equation governs the steady-state pressure distribution,
\[{{\partial^{2} p} \over {\partial x^{2}}} + {{\partial^{2} p} \over {\partial y^{2}}} + {{\partial^{2} p} \over {\partial z^{2}}} = 0\]
The solution is given by,
\[p(x,y,z) = p_{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.¶
Comparison of the PFLOTRAN vs. analytical solution for RICHARDS 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 10
DXYZ
0.1d0
0.1d0
0.1d0
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.5
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-10
PERM_Y 1.d-10
PERM_Z 1.d-10
/
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.05d0
M 0.5d0
/
PERMEABILITY_FUNCTION MUALEM_VG_GAS
LIQUID_RESIDUAL_SATURATION 0.05d0
GAS_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
END
STRATA
REGION all
MATERIAL cube
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.10 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET pressure_bc_west
HDF5_DATASET_NAME node_centered_surf_west
FILENAME ../dataset.h5
END
DATASET pressure_bc_east
HDF5_DATASET_NAME node_centered_surf_east
FILENAME ../dataset.h5
END
DATASET pressure_bc_north
HDF5_DATASET_NAME node_centered_surf_north
FILENAME ../dataset.h5
END
DATASET pressure_bc_south
HDF5_DATASET_NAME node_centered_surf_south
FILENAME ../dataset.h5
END
DATASET pressure_bc_bottom
HDF5_DATASET_NAME node_centered_surf_bottom
FILENAME ../dataset.h5
END
DATASET pressure_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 25.d0 C
LIQUID_PRESSURE 1.0 MPa
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_west
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_east
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_north
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_south
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION bottom_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_bottom
MOLE_FRACTION 1.d-10
END
FLOW_CONDITION top_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
MOLE_FRACTION dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_top
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 10 10 10
DXYZ
0.1d0
0.1d0
0.1d0
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.5
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-10
PERM_Y 1.d-10
PERM_Z 1.d-10
/
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
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.10 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET pressure_bc_west
HDF5_DATASET_NAME node_centered_surf_west
FILENAME ../dataset.h5
END
DATASET pressure_bc_east
HDF5_DATASET_NAME node_centered_surf_east
FILENAME ../dataset.h5
END
DATASET pressure_bc_north
HDF5_DATASET_NAME node_centered_surf_north
FILENAME ../dataset.h5
END
DATASET pressure_bc_south
HDF5_DATASET_NAME node_centered_surf_south
FILENAME ../dataset.h5
END
DATASET pressure_bc_bottom
HDF5_DATASET_NAME node_centered_surf_bottom
FILENAME ../dataset.h5
END
DATASET pressure_bc_top
HDF5_DATASET_NAME node_centered_surf_top
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE 1.0 MPa
END
FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_west
END
FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_east
END
FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_north
END
FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_south
END
FLOW_CONDITION bottom_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_bottom
END
FLOW_CONDITION top_face
TYPE
TEMPERATURE dirichlet
LIQUID_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
LIQUID_PRESSURE DATASET pressure_bc_top
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 (RICHARDS Mode)¶
The RICHARDS Mode PFLOTRAN input file can be downloaded
here.
SIMULATION
SIMULATION_TYPE SUBSURFACE
PROCESS_MODELS
SUBSURFACE_FLOW flow
MODE RICHARDS
/
/
END
SUBSURFACE
GRID
TYPE structured
NXYZ 10 10 10
DXYZ
0.1d0
0.1d0
0.1d0
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.5
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
CHARACTERISTIC_CURVES cc1
PERMEABILITY
PERM_X 1.d-10
PERM_Y 1.d-10
PERM_Z 1.d-10
/
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.05d0
M 0.5d0
/
PERMEABILITY_FUNCTION MUALEM_VG_GAS
LIQUID_RESIDUAL_SATURATION 0.05d0
GAS_RESIDUAL_SATURATION 0.1d0
M 0.5d0
/
END
STRATA
REGION all
MATERIAL cube
END
FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END
EOS WATER
DENSITY CONSTANT 1000.d0 kg/m^3
VISCOSITY CONSTANT 1.0d-3 Pa-s
END
TIME
FINAL_TIME 1 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.10 yr at 0.d0 yr
END
OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
FORMAT HDF5
/
END
DATASET pressure_bc_west
HDF5_DATASET_NAME node_centered_surf_west
FILENAME ../dataset.h5
END
DATASET pressure_bc_east
HDF5_DATASET_NAME node_centered_surf_east
FILENAME ../dataset.h5
END
DATASET pressure_bc_north
HDF5_DATASET_NAME node_centered_surf_north
FILENAME ../dataset.h5
END
DATASET pressure_bc_south
HDF5_DATASET_NAME node_centered_surf_south
FILENAME ../dataset.h5
END
DATASET pressure_bc_bottom
HDF5_DATASET_NAME node_centered_surf_bottom
FILENAME ../dataset.h5
END
DATASET pressure_bc_top
HDF5_DATASET_NAME node_centered_surf_top
FILENAME ../dataset.h5
END
FLOW_CONDITION initial
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE 1.0 MPa
END
FLOW_CONDITION west_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_west
END
FLOW_CONDITION east_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_east
END
FLOW_CONDITION north_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_north
END
FLOW_CONDITION south_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_south
END
FLOW_CONDITION bottom_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_bottom
END
FLOW_CONDITION top_face
TYPE
LIQUID_PRESSURE dirichlet
/
LIQUID_PRESSURE DATASET pressure_bc_top
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')
p0 = 1.e6 # [Pa]
L = 1. # [m]
# 2D Surface:
# -----------------------------------------------------------------------------
# Pressure 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] = p0*(0+(y/L)+(z/L))
h5dset = h5grp.create_dataset('Data', data=rarray)
# 2D Surface:
# -----------------------------------------------------------------------------
# Pressure 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] = p0*(1+(y/L)+(z/L))
h5dset = h5grp.create_dataset('Data', data=rarray)
# 2D Surface:
# -----------------------------------------------------------------------------
# Pressure 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] = p0*((x/L)+(y/L)+0)
h5dset = h5grp.create_dataset('Data', data=rarray)
# 2D Surface:
# -----------------------------------------------------------------------------
# Pressure 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] = p0*((x/L)+(y/L)+1)
h5dset = h5grp.create_dataset('Data', data=rarray)
# 2D Surface:
# -----------------------------------------------------------------------------
# Pressure 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] = p0*((x/L)+0+(z/L))
h5dset = h5grp.create_dataset('Data', data=rarray)
# 2D Surface:
# -----------------------------------------------------------------------------
# Pressure 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] = p0*((x/L)+1+(z/L))
h5dset = h5grp.create_dataset('Data', data=rarray)
The Python Script¶
def flow_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.2.5, pg.31
# "A 3D Steady-State Pressure Distribution"
#
# Author: Jenn Frederick
# Date: 07/25/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] = 1. # [m] lx
lxyz[1] = 1. # [m] ly
lxyz[2] = 1. # [m] lz
nxyz[0] = 10 # [m] nx
nxyz[1] = 10 # [m] ny
nxyz[2] = 10 # [m] nz
dxyz = lxyz/nxyz # [m]
p0 = 1. # [MPa]
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
p_soln = np.zeros((nx,ny,nz)) # [C]
p_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
p_soln[i,j,k] = p0*((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+00 y/Liquid_Pressure [Pa]'
p_pflotran[:,:,:] = read_pflotran_output_3D(path+
'/3D_steady_pressure_BC_1st_kind.h5',index_string,remove)
ierr = check(p_pflotran)
# Convert units [Pa -> MPa]
p_pflotran = p_pflotran/1.e6
max_percent_error = calc_relative_error(p_soln,p_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 = nxyz*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,p_soln,p_pflotran,'Distance [m]',
'Pressure [MPa]',"{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.
