1D Steady Gas (Pressure), BCs of 1st Kind¶

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The Python Script

The Problem Description¶

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

The domain is a 40x5x5 meter rectangular beam extending along the positive x-axis and is made up of 40x1x1 cubic grid cells with dimensions 1x5x5 meters. The materials are assigned the following properties: porosity = 0.5; gas viscosity $$\mu$$ = 1.0e-5 Pa-s; permeability = 1e-15 m^2; gas saturation $$S_g$$ = 1.0 in the entire domain.

The gas pressure is initially uniform at p(t=0) = 1.5e5 Pa. The boundary gas pressures are p(x=0) = 2.0e5 Pa and p(x=L) = 1.0e5 Pa, where L = 40 m. The simulation is run until the steady-state pressure distribution develops.

The LaPlace equation governs the steady-state gas pressure distribution,

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

The solution is given by,

$p(x) = \sqrt{ax + b}$

When the gas pressure boundary conditions are applied, the solution becomes,

$p(x) = \sqrt{\left({p_1^2 - p_0^2}\right) {x \over L} + p_0^2}$

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
4.0d0
1.0d0
1.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
40.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
40.0d0 0.d0 0.d0
40.0d0 1.0d0 1.0d0
/
END

MATERIAL_PROPERTY beam
ID 1
POROSITY 0.50
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 100.
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.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 beam
END

EOS GAS
VISCOSITY CONSTANT 1.0d-5 Pa-s
#DENSITY CONSTANT 1.0 kg/m^3
END

FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
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
GAS_SATURATION dirichlet
GAS_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
GAS_SATURATION 1.0
GAS_PRESSURE 1.5d5 Pa
END

FLOW_CONDITION left_end
TYPE
TEMPERATURE dirichlet
GAS_SATURATION dirichlet
GAS_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
GAS_SATURATION 1.0
GAS_PRESSURE 2.0d5 Pa
END

FLOW_CONDITION right_end
TYPE
TEMPERATURE dirichlet
GAS_SATURATION dirichlet
GAS_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
GAS_SATURATION 1.0
GAS_PRESSURE 1.0d5 Pa
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
REGION left_end
FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
FLOW_CONDITION right_end
REGION right_end
END

END_SUBSURFACE


The Python Script¶

def gas_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.3.1, pg.40
# "A 1D Steady-State Gas Pressure Distribution, Boundary Conditions of 1st Kind"
#
# Author: Jenn Frederick
# Date: 09/09/2016
#*******************************************************************************
nxyz = np.zeros(3) + 1
dxyz = np.zeros(3) + 1.
lxyz = np.zeros(3) + 1.
error_analysis = np.zeros(num_tries)
dxyz_record = np.zeros((num_tries,3))
test_pass = False
try_count = 0

# initial discretization values
lxyz[0] = 40.     # [m] lx
nxyz[0] = 10      # [m] nx
dxyz = lxyz/nxyz  # [m]
p0 = 2.0e5        # [Pa]
p1 = 1.0e5        # [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]
x_soln = np.concatenate((x_soln,[Lx]),axis=0)   # [m]
p_soln = np.zeros(nx+2)                         # [Pa]
p_pflotran = np.zeros(nx)                       # [Pa]

# create the analytical solution
p_soln = np.sqrt((((p1*p1) - (p0*p0))*(x_soln/Lx)) + (p0*p0)) # [Pa]

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

index_string = 'Time:  1.00000E+00 y/Gas_Pressure [Pa]'
ierr = check(p_pflotran)

# convert pressure units [Pa -> MPa]
p_pflotran = p_pflotran/1.0e6         # [MPa]
p_soln = p_soln/1.0e6                 # [MPa]

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

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

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

# Add test result to report card

return;


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

1D Steady Gas (Pressure), BCs of 1st & 2nd Kind¶

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The Python Script

The Problem Description¶

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

The domain is a 40x5x5 meter rectangular beam extending along the positive x-axis and is made up of 40x1x1 cubic grid cells with dimensions 1x5x5 meters. The materials are assigned the following properties: porosity = 0.5; gas viscosity $$\mu$$ = 1.0e-5 Pa-s; permeability = 1e-15 m^2; gas saturation $$S_g$$ = 1.0 in the entire domain.

The gas pressure is initially uniform at p(t=0) = 1.0e5 Pa. At the east boundary, a dirichlet gas pressure p(x=L) = 1.0e5 Pa is applied, where L = 40 m. At the west boundary, a specific gas flow is applied Q(x=0) = 0.17 Pa m/s,

$Q = -\frac{k}{\mu}p\frac{dp}{dx}$

PFLOTRAN does not have the capability to handle such a boundary condition. To convert specific gas flow to a gas flux in m/s, the ideal gas law is used by dividing the specific gas flow by the boundary pressure,

$p = \frac{n}{V}RT$

where $$\frac{n}{V}$$ is the gas density at the boundary in mol m-1, R is the ideal gas constant (8.314), and T is the temperature in Kelvin. Therefore, the specific gas flow boundary at x=0 is converted to a gas flux boundary as,

(1)$q = \frac{0.17}{\frac{n}{V}RT}$

As gas flows into the domain, pressure increases, and as a result gas density also increases. To figure out what value of gas flux must be applied at the west boundary, we can use the analytical solution for gas pressure at x=0. Assuming a certain temperature, the gas density can be calculated using the ideal gas law. That gas density can then be used in equation (1) to set the value of the gas flux boundary condition. In this example, the gas pressure at steady state at the west boundary is approximately 4.4358 kg m-3, which is converted to approximately 153.17 mol m-3, giving $$\frac{n}{V}$$. The values for R and T are known.

The simulation is run until the steady-state pressure distribution develops, which is 20 years.

The LaPlace equation governs the steady-state gas pressure distribution,

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

The solution is given by,

$p(x) = \sqrt{ax + b}$

When the boundary conditions are applied, the solution becomes,

$p(x) = \sqrt{\frac{2 Q \mu}{k}(L-x) + p_1^2}$

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
2.0d0
1.0d0
1.0d0
END
END

REGION all
COORDINATES
0.0d0 0.0d0 0.0d0
40.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
40.0d0 0.d0 0.d0
40.0d0 1.0d0 1.0d0
/
END

REGION left_grid_cell
COORDINATE 1.0d0 0.5d0 0.5d0
END

MATERIAL_PROPERTY beam
ID 1
POROSITY 0.50
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 100.
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.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 beam
END

EOS GAS
VISCOSITY CONSTANT 1.0d-5 Pa-s
END

FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END

TIME
FINAL_TIME 20 yr
INITIAL_TIMESTEP_SIZE 1.d-4 day
MAXIMUM_TIMESTEP_SIZE 0.01d0 yr at 0.d0 yr
END

OBSERVATION
REGION left_grid_cell
END

OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
MOLE_FRACTION dirichlet
GAS_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
MOLE_FRACTION 0.98
GAS_PRESSURE 1.0d5 Pa
END

FLOW_CONDITION left_end
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0.d0 m/s
# 0.17/((n/V)*R*T) = 0.17/(153.17*8.314*298.15)
GAS_FLUX 4.475d-7 m/s
END

FLOW_CONDITION right_end
TYPE
TEMPERATURE dirichlet
MOLE_FRACTION dirichlet
GAS_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
MOLE_FRACTION 0.98
GAS_PRESSURE 1.0d5 Pa
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION left_end
REGION left_end
FLOW_CONDITION left_end
END

BOUNDARY_CONDITION right_end
FLOW_CONDITION right_end
REGION right_end
END

END_SUBSURFACE


The Python Script¶

def gas_steady_1D_BC1st2ndkind(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.3.2, pg.41
# "A 1D Steady-State Gas Pressure Distribution, Boundary Conditions of 1st
# and 2nd Kind"
#
# Author: Jenn Frederick
# Date: 11/03/2016
#*******************************************************************************
nxyz = np.zeros(3) + 1
dxyz = np.zeros(3) + 1.
lxyz = np.zeros(3) + 1.
error_analysis = np.zeros(num_tries)
dxyz_record = np.zeros((num_tries,3))
test_pass = False
try_count = 0

# initial discretization values
lxyz[0] = 40.     # [m] lx
nxyz[0] = 20      # [m] nx
dxyz = lxyz/nxyz  # [m]
k = 1.0e-15    # [m^2]
mu = 1.0e-5    # [Pa*s]
p1 = 1.0e5     # [Pa]
Q = 0.17       # [Pa*m/s]
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)                           # [Pa]
p_pflotran = np.zeros(nx)                         # [Pa]

# create the analytical solution
p_soln = np.sqrt(((2*Q*mu)/k)*(Lx-x_soln)+(p1**2)) # [Pa]

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

index_string = 'Time:  2.00000E+01 y/Gas_Pressure [Pa]'
ierr = check(p_pflotran)

# convert pressure units [Pa -> MPa]
p_pflotran = p_pflotran/1.0e6         # [MPa]
p_soln = p_soln/1.0e6                 # [MPa]

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

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

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

# Add test result to report card

return;


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

2D Steady Gas (Pressure), BCs of 1st & 2nd Kind¶

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The Dataset

The Python Script

The Problem Description¶

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

The domain is a 1x1x0.5 meter rectangular plate extending along the positive x-axis and y-axis and is made up of 20x20x1 hexahedral grid cells with dimensions 0.05x0.05x0.5 meters. The materials are assigned the following properties: porosity = 0.5; gas viscosity $$\mu$$ = 1.0e-5 Pa-s; permeability = 1e-15 m^2; gas saturation $$S_g$$ = 1.0 in the entire domain.

The gas pressure is initially uniform at p(t=0) = 1.0e5 Pa. At the west and south boundaries, a dirichlet gas pressure p(x=0,y=(0,L)) = p(x=(0,L),y=0) = 1.0e5 Pa is applied, where L = 1 m. At the north and east boundaries, a specific gas flow [Pa m/s] is applied that is a function of x and y,

\begin{align}\begin{aligned}Q = \frac{3 k p_0^2 y}{2 \mu L^2} \hspace{0.5 in} east \hspace{0.2 in} face\\Q = \frac{3 k p_0^2 x}{2 \mu L^2} \hspace{0.5 in} north \hspace{0.2 in} face\end{aligned}\end{align}

PFLOTRAN does not have the capability to handle such a boundary condition. For an explanation on how a specific gas flow is treated, see 1D Steady Gas (Pressure), BCs of 1st & 2nd Kind.

The simulation is run until the steady-state pressure distribution develops, which is 20 years.

The LaPlace equation governs the steady-state gas pressure distribution,

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

When the boundary conditions are applied, the solution becomes,

$p(x,y) = p_0 \sqrt{1 + 3\frac{xy}{L^2}}$

where $$p_0$$ is the initial pressure, p(t=0) = 1.0e5 Pa.

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 15 15 1
DXYZ
0.0666666666667d0
0.0666666666667d0
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.50
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 100.
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.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 plate
END

FLUID_PROPERTY
DIFFUSION_COEFFICIENT 1.d-9
END

EOS GAS
VISCOSITY CONSTANT 1.0d-5 Pa-s
END

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

OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END

DATASET flux_bc_north
HDF5_DATASET_NAME x_line_cell_centered_north
FILENAME ../dataset.h5
END

DATASET flux_bc_east
HDF5_DATASET_NAME y_line_cell_centered_east
FILENAME ../dataset.h5
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
MOLE_FRACTION dirichlet
GAS_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
MOLE_FRACTION 0.98
GAS_PRESSURE 1.0d5 Pa
END

FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0.d0 m/s
GAS_FLUX DATASET flux_bc_east
END

FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0.d0 m/s
GAS_FLUX DATASET flux_bc_north
END

FLOW_CONDITION dirichlet_pressure
TYPE
TEMPERATURE dirichlet
MOLE_FRACTION dirichlet
GAS_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
MOLE_FRACTION 0.98
GAS_PRESSURE 1.0d5 Pa
END

INITIAL_CONDITION initial
REGION all
FLOW_CONDITION initial
END

BOUNDARY_CONDITION west_face
REGION west_face
FLOW_CONDITION dirichlet_pressure
END

BOUNDARY_CONDITION south_face
REGION south_face
FLOW_CONDITION dirichlet_pressure
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 boundary conditions is created with the following python script called create_dataset.py:

import sys
from h5py import *
import numpy as np

nx = 1; dx = 1; lx = 1
ny = 1; dy = 1; ly = 1
nz = 1; dz = 1; lz = 1
for k in range(len(sys.argv)-1):
k = k + 1
if sys.argv[k] == '-nx':
nx = int(float(sys.argv[k+1]))
if sys.argv[k] == '-ny':
ny = int(float(sys.argv[k+1]))
if sys.argv[k] == '-nz':
nz = int(float(sys.argv[k+1]))
if sys.argv[k] == '-dx':
dx = float(sys.argv[k+1])
if sys.argv[k] == '-dy':
dy = float(sys.argv[k+1])
if sys.argv[k] == '-dz':
dz = float(sys.argv[k+1])
if sys.argv[k] == '-lx':
lx = float(sys.argv[k+1])
if sys.argv[k] == '-ly':
ly = float(sys.argv[k+1])
if sys.argv[k] == '-lz':
lz = float(sys.argv[k+1])

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

p0 = 1.0e5   # [Pa]
k = 1.0e-15  # [m2]
mu = 1.0e-5  # [Pa-s]

# 1d line in x
# Flux boundary condition; NORTH
h5grp = h5file.create_group('x_line_cell_centered_north')
h5grp.attrs['Cell Centered'] = True
h5grp.attrs['Interpolation Method'] = np.string_('STEP')
h5grp.attrs['Dimension'] = np.string_('X')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dx]
# Location of origin
h5grp.attrs['Origin'] = [0.]
rarray = np.zeros(nx,'=f8')
for i in range(nx):
x = (float(i)*lx)/nx + dx/2.
p_gas = p0*np.sqrt(1+3*(x/lx))  # [Pa] from analytical solution at y=L
rarray[i] = ((3*k)/(2*mu))*((p0**2)/lx)*(x/lx)   # [Pa m/s]
rarray[i] = rarray[i]/p_gas                      # [m/s]
h5dset = h5grp.create_dataset('Data', data=rarray)

# 1d line in y
# Flux boundary condition; EAST
h5grp = h5file.create_group('y_line_cell_centered_east')
h5grp.attrs['Cell Centered'] = True
h5grp.attrs['Interpolation Method'] = np.string_('STEP')
h5grp.attrs['Dimension'] = np.string_('Y')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dy]
# Location of origin
h5grp.attrs['Origin'] = [0.]
rarray = np.zeros(ny,'=f8')
for i in range(ny):
y = (float(i)*ly)/ny + dy/2.
p_gas = p0*np.sqrt(1+3*(y/ly))  # [Pa] from analytical solution at x=L
rarray[i] = ((3*k)/(2*mu))*((p0**2)/ly)*(y/ly)   # [Pa m/s]
rarray[i] = rarray[i]/p_gas                      # [m/s]
h5dset = h5grp.create_dataset('Data', data=rarray)


The Python Script¶

def gas_steady_2D_BC1st2ndkind(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.3.3, pg.42
# "A 2D Steady-State Gas Pressure Distribution, Boundary Conditions of 1st
# and 2nd Kind"
#
# Author: Jenn Frederick
# Date: 11/04/2016
#*******************************************************************************
nxyz = np.zeros(3) + 1
dxyz = np.zeros(3) + 1.
lxyz = np.zeros(3) + 1.
error_analysis = np.zeros(num_tries)
dxyz_record = np.zeros((num_tries,3))
test_pass = False
try_count = 0

# initial discretization values
lxyz[0] = 1.      # [m] lx
lxyz[1] = 1.      # [m] ly
nxyz[0] = 15      # [m] nx
nxyz[1] = 15      # [m] ny
dxyz = lxyz/nxyz  # [m]
p0 = 1.0e5        # [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

p_soln = np.zeros((nx,ny))                       # [Pa]
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*np.sqrt(1+3*((x*y)/(Lx*Ly)))

# run PFLOTRAN simulation
os.chdir('../')
command_line = 'python create_dataset.py'
command_line = command_line + ' -nx ' + str(nx) + ' -ny ' + str(ny)
command_line = command_line + ' -lx ' + str(Lx) + ' -ly ' + str(Ly)
command_line = command_line + ' -dx ' + str(dx) + ' -dy ' + str(dy)
os.system(command_line)
os.chdir('general_mode')
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)

index_string = 'Time:  1.00000E+01 y/Gas_Pressure [Pa]'
ierr = check(p_pflotran)

# convert pressure units [Pa -> MPa]
p_pflotran = p_pflotran/1.0e6         # [MPa]
p_soln = p_soln/1.0e6                 # [MPa]

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

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

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

# Add test result to report card

return;


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

3D Steady Gas (Pressure), BCs of 2nd Kind¶

The Problem Description

The PFLOTRAN Input File (GENERAL Mode)

The Dataset

The Python Script

The Problem Description¶

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

The domain is a 1x1x1 meter cube extending along the positive x-axis, y-axis, and z-axis and is made up of 20x20x20 hexahedral grid cells with dimensions 0.05x0.05x0.05 meters. The material is assigned the following properties: porosity = 0.5; gas viscosity $$\mu$$ = 1.0e-5 Pa-s; permeability = 1e-15 m^2; gas saturation $$S_g$$ = 1.0 in the entire domain.

The gas pressure is initially uniform at p(t=0) = 1.0e5 Pa. At each face of the cube, a specific gas flow [Pa m/s] boundary condition is applied:

\begin{align}\begin{aligned}-Q(0,y,z) = Q(L,y,z) = \frac{3 k p_0^2 y}{4 \mu L^2}\\-Q(x,0,z) = Q(x,L,z) = \frac{3 k p_0^2 x}{4 \mu L^2}\\-Q(x,y,0) = Q(x,y,L) = \frac{3 k p_0^2 }{4 \mu L}\end{aligned}\end{align}

Note that this specific gas flow boundary condition acts as a source on the faces where x, y, or z = L, and a sink on the faces where x, y, or z = 0. PFLOTRAN does not have the capability to handle such a boundary condition. For an explanation on how a specific gas flow is treated, see 1D Steady Gas (Pressure), BCs of 1st & 2nd Kind.

The simulation is run until the steady-state pressure distribution develops, which is 10 years.

The LaPlace equation governs the steady-state gas pressure distribution,

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

When the boundary conditions are applied, the solution becomes,

$p(x,y,z) = p_0 \sqrt{1 + \frac{3}{2}\left({\frac{xy}{L^2}+\frac{z}{L}}\right)}$

where $$p_0$$ is the initial pressure, p(t=0) = 1.0e5 Pa.

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 20 20
DXYZ
0.05d0
0.05d0
0.05d0
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.50
TORTUOSITY 1.d0
ROCK_DENSITY 2.8E3
HEAT_CAPACITY 100.
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.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 GAS
VISCOSITY CONSTANT 1.0d-5 Pa-s
END

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

OUTPUT
SNAPSHOT_FILE
NO_PRINT_INITIAL
#NO_PRINT_FINAL
FORMAT HDF5
/
END

DATASET flux_bc_west
HDF5_DATASET_NAME cell_centered_surf_west
FILENAME ../dataset.h5
END

DATASET flux_bc_east
HDF5_DATASET_NAME cell_centered_surf_east
FILENAME ../dataset.h5
END

DATASET flux_bc_north
HDF5_DATASET_NAME cell_centered_surf_north
FILENAME ../dataset.h5
END

DATASET flux_bc_south
HDF5_DATASET_NAME cell_centered_surf_south
FILENAME ../dataset.h5
END

DATASET flux_bc_top
HDF5_DATASET_NAME cell_centered_surf_top
FILENAME ../dataset.h5
END

DATASET flux_bc_bottom
HDF5_DATASET_NAME cell_centered_surf_bottom
FILENAME ../dataset.h5
END

FLOW_CONDITION initial
TYPE
TEMPERATURE dirichlet
MOLE_FRACTION dirichlet
GAS_PRESSURE dirichlet
/
TEMPERATURE 25.d0 C
MOLE_FRACTION 0.98
GAS_PRESSURE 1.0d5 Pa
END

FLOW_CONDITION east_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0.d0 m/s
GAS_FLUX DATASET flux_bc_east
END

FLOW_CONDITION west_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0.d0 m/s
GAS_FLUX DATASET flux_bc_west
END

FLOW_CONDITION north_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0.d0 m/s
GAS_FLUX DATASET flux_bc_north
END

FLOW_CONDITION south_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0.d0 m/s
GAS_FLUX DATASET flux_bc_south
END

FLOW_CONDITION top_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0.d0 m/s
GAS_FLUX DATASET flux_bc_top
END

FLOW_CONDITION bottom_face
TYPE
TEMPERATURE dirichlet
LIQUID_FLUX neumann
GAS_FLUX neumann
/
TEMPERATURE 25.d0 C
LIQUID_FLUX 0.d0 m/s
GAS_FLUX DATASET flux_bc_bottom
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 boundary conditions is created with the following python script called create_dataset.py:

import sys
from h5py import *
import numpy as np

nx = 1; dx = 1; lx = 1
ny = 1; dy = 1; ly = 1
nz = 1; dz = 1; lz = 1
for k in range(len(sys.argv)-1):
k = k + 1
if sys.argv[k] == '-nx':
nx = int(float(sys.argv[k+1]))
if sys.argv[k] == '-ny':
ny = int(float(sys.argv[k+1]))
if sys.argv[k] == '-nz':
nz = int(float(sys.argv[k+1]))
if sys.argv[k] == '-dx':
dx = float(sys.argv[k+1])
if sys.argv[k] == '-dy':
dy = float(sys.argv[k+1])
if sys.argv[k] == '-dz':
dz = float(sys.argv[k+1])
if sys.argv[k] == '-lx':
lx = float(sys.argv[k+1])
if sys.argv[k] == '-ly':
ly = float(sys.argv[k+1])
if sys.argv[k] == '-lz':
lz = float(sys.argv[k+1])

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

p0 = 1.0e5   # [Pa]
K = 1.0e-15  # [m2]
mu = 1.0e-5  # [Pa-s]

# 2D Surface:
# -----------------------------------------------------------------------------
# Gas flux boundary condition x=0 face
h5grp = h5file.create_group('cell_centered_surf_west')
h5grp.attrs['Dimension'] = np.string_('YZ')
h5grp.attrs['Cell Centered'] = True
h5grp.attrs['Interpolation Method'] = np.string_('STEP')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dy,dz]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
rarray = np.zeros((ny,nz),'=f8')
for j in range(ny):
for k in range(nz):
x = 0.
y = (float(j)*ly)/(ny) + dy/2.
z = (float(k)*lz)/(nz) + dz/2.
p_gas = p0*np.sqrt(1.+(3./2.)*((x/lx)*(y/ly)+(z/lz)))
rarray[j][k] = ((3.*K)/(4.*mu))*(p0**2/lx)*(y/ly)
# flow is out: (-)
rarray[j][k] = -1.0*rarray[j][k]/p_gas
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Gas flux boundary condition x=L face
h5grp = h5file.create_group('cell_centered_surf_east')
h5grp.attrs['Dimension'] = np.string_('YZ')
h5grp.attrs['Cell Centered'] = True
h5grp.attrs['Interpolation Method'] = np.string_('STEP')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dy,dz]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
rarray = np.zeros((ny,nz),'=f8')
for j in range(ny):
for k in range(nz):
x = lx
y = (float(j)*ly)/(ny) + dy/2.
z = (float(k)*lz)/(nz) + dz/2.
p_gas = p0*np.sqrt(1+(3/2)*((x/lx)*(y/ly)+(z/lz)))  # [Pa] from analytical solution
rarray[j][k] = ((3.*K)/(4.*mu))*(p0**2/lx)*(y/ly)   # [Pa m/s]
rarray[j][k] = rarray[j][k]/p_gas                   # [m/s]
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Gas flux boundary condition z=0 face
h5grp = h5file.create_group('cell_centered_surf_bottom')
h5grp.attrs['Dimension'] = np.string_('XY')
h5grp.attrs['Cell Centered'] = True
h5grp.attrs['Interpolation Method'] = np.string_('STEP')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dx,dy]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
rarray = np.zeros((nx,ny),'=f8')
for i in range(nx):
for j in range(ny):
x = (float(i)*lx)/(nx) + dx/2.
y = (float(j)*ly)/(ny) + dy/2.
z = 0.
p_gas = p0*np.sqrt(1.+(3./2.)*((x/lx)*(y/ly)+(z/lz)))
rarray[i][j] = ((3.*K)/(4.*mu))*(p0**2/lz)
# flow is out: (-)
rarray[i][j] = -1.0*rarray[i][j]/p_gas
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Gas flux boundary condition z=L face
h5grp = h5file.create_group('cell_centered_surf_top')
h5grp.attrs['Dimension'] = np.string_('XY')
h5grp.attrs['Cell Centered'] = True
h5grp.attrs['Interpolation Method'] = np.string_('STEP')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dx,dy]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
rarray = np.zeros((nx,ny),'=f8')
for i in range(nx):
for j in range(ny):
x = (float(i)*lx)/(nx) + dx/2.
y = (float(j)*ly)/(ny) + dy/2.
z = lz
p_gas = p0*np.sqrt(1.+(3./2.)*((x/lx)*(y/ly)+(z/lz)))
rarray[i][j] = ((3.*K)/(4.*mu))*(p0**2/lz)
rarray[i][j] = rarray[i][j]/p_gas
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Gas flux boundary condition y=0 face
h5grp = h5file.create_group('cell_centered_surf_south')
h5grp.attrs['Dimension'] = np.string_('XZ')
h5grp.attrs['Cell Centered'] = True
h5grp.attrs['Interpolation Method'] = np.string_('STEP')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dx,dz]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
rarray = np.zeros((nx,nz),'=f8')
for i in range(nx):
for k in range(nz):
x = (float(i)*lx)/(nx) + dx/2.
y = 0.
z = (float(k)*lz)/(nz) + dz/2.
p_gas = p0*np.sqrt(1.+(3./2.)*((x/lx)*(y/ly)+(z/lz)))
rarray[i][k] = x*((3.*K)/(4.*mu))*(p0**2/ly**2)
# flow is out: (-)
rarray[i][k] = -1.0*rarray[i][k]/p_gas
h5dset = h5grp.create_dataset('Data', data=rarray)

# 2D Surface:
# -----------------------------------------------------------------------------
# Gas flux boundary condition y=L face
h5grp = h5file.create_group('cell_centered_surf_north')
h5grp.attrs['Dimension'] = np.string_('XZ')
h5grp.attrs['Cell Centered'] = True
h5grp.attrs['Interpolation Method'] = np.string_('STEP')
# Delta length between points [m]
h5grp.attrs['Discretization'] = [dx,dz]
# Location of origin
h5grp.attrs['Origin'] = [0.,0.]
rarray = np.zeros((nx,nz),'=f8')
for i in range(nx):
for k in range(nz):
x = (float(i)*lx)/(nx) + dx/2.
y = ly
z = (float(k)*lz)/(nz) + dz/2.
p_gas = p0*np.sqrt(1.+(3./2.)*((x/lx)*(y/ly)+(z/lz)))
rarray[i][k] = x*((3.*K)/(4.*mu))*(p0**2/ly**2)
rarray[i][k] = rarray[i][k]/p_gas
h5dset = h5grp.create_dataset('Data', data=rarray)


The Python Script¶

def gas_steady_3D_BC2ndkind(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.3.4, pg.43
# "A 3D Steady-State Gas Pressure Distribution, Boundary Conditions of 2nd
# Kind"
#
# Author: Jenn Frederick
# Date: 11/07/2016
#*******************************************************************************
nxyz = np.zeros(3) + 1
dxyz = np.zeros(3) + 1.
lxyz = np.zeros(3) + 1.
error_analysis = np.zeros(num_tries)
dxyz_record = np.zeros((num_tries,3))
test_pass = False
try_count = 0

# initial discretization values
lxyz[0] = 1.      # [m] lx
lxyz[1] = 1.      # [m] ly
lxyz[2] = 1.      # [m] lz
nxyz[0] = 5      # [m] nx
nxyz[1] = 5      # [m] ny
nxyz[2] = 5      # [m] nz
dxyz = lxyz/nxyz  # [m]
p0 = 1.0e5     # [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

p_soln = np.zeros((nx,ny,nz))                    # [Pa]
p_pflotran = np.zeros((nx,ny,nz))                # [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]
z_soln = np.linspace(0.+(dz/2.),Lz-(dz/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*np.sqrt(1.+(3./2.)*((x/Lx)*(y/Ly)+(z/Lz)))

# run PFLOTRAN simulation
os.chdir('../')
command_line = 'python create_dataset.py'
command_line = command_line + ' -nx ' +str(nx)+ ' -ny ' +str(ny)+ ' -nz ' +str(nz)
command_line = command_line + ' -lx ' +str(Lx)+ ' -ly ' +str(Ly)+ ' -lz ' +str(Lz)
command_line = command_line + ' -dx ' +str(dx)+ ' -dy ' +str(dy)+ ' -dz ' +str(dz)
os.system(command_line)
os.chdir('general_mode')
run_pflotran(input_prefix,nxyz,dxyz,lxyz,remove,screen_on,pf_exe,mpi_exe)

index_string = 'Time:  1.00000E+01 y/Gas_Pressure [Pa]'
ierr = check(p_pflotran)

# convert pressure units [Pa -> MPa]
p_pflotran = p_pflotran/1.0e6         # [MPa]
p_soln = p_soln/1.0e6                 # [MPa]

max_percent_error = calc_relative_error(p_soln,p_pflotran,ierr)
record_error(error_analysis,dxyz_record,max_percent_error,dxyz,try_count)
test_pass = does_pass(max_percent_error,try_count,num_tries)
nxyz = nxyz*2.
dxyz = lxyz/nxyz

# Plot the PFLOTRAN and analytical solutions

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.