Porous Flow Tutorial Page 06. Adding a tracer
On this Page we will depart from PorousFlowBasicTHM and use PorousFlowFullySaturated. This action employs mass lumping but no numerical stabilization by default. Therefore, the results will be slightly different than those achieved by PorousFlowBasicTHM (no lumping or stabilization) and from the remainder of PorousFlow (mass lumping and either full upwinding or TVD stabilization).
The reason for using PorousFlowFullySaturated is that it allows multi-component physics to be studied. Denoting the mass fraction by (where ), the fluid equations are (1) In this equation the fluid density, , is bracketed, because the user has the option of including it or not using the multiply_by_density flag. The advantage of setting multiply_by_density = false is that the equations become more linear, but it means that care must be taken when using other parts of PorousFlow. For instance, the PorousFlowMass Postprocessor is measuring fluid mass, not fluid volume.
Tracers may now be incorporated into the simulation. A tracer_concentration variable is defined, and set to initial conditions where it is 0.5 at the injection area, and zero elsewhere. Modifying the isothermal simulation of Page 01, PorousFlowFullySaturated is provided with its name:
[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
[porepressure]
[]
[tracer_concentration]
[]
[]
[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
[tracer_concentration]
type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
function<<<{"description": "The initial condition function."}>>> = '0.5*if(x*x+y*y<1.01,1,0)'
variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = tracer_concentration
[]
[]
[PorousFlowFullySaturated<<<{"href": "../../syntax/PorousFlowFullySaturated/index.html"}>>>]
porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
coupling_type<<<{"description": "The type of simulation. For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient. For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material"}>>> = Hydro
gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = the_simple_fluid
mass_fraction_vars<<<{"description": "List of variables that represent the mass fractions. With only one fluid component, this may be left empty. With N fluid components, the format is 'f_0 f_1 f_2 ... f_(N-1)'. That is, the N^th component need not be specified because f_N = 1 - (f_0 + f_1 + ... + f_(N-1)). It is best numerically to choose the N-1 mass fraction variables so that they represent the fluid components with small concentrations. This Action will associated the i^th mass fraction variable to the equation for the i^th fluid component, and the pressure variable to the N^th fluid component."}>>> = tracer_concentration
stabilization<<<{"description": "Numerical stabilization used. 'Full' means full upwinding. 'KT' means FEM-TVD stabilization of Kuzmin-Turek"}>>> = none # Note to reader: 06_KT.i uses KT stabilization - compare the results
[]To encourage flow, rather than a rapid convergence to a more-or-less steady-state situation, the outer boundary is fixed to . The tracer concentration is fixed to 0.5 at the injection area.
[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
[constant_injection_porepressure]
type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = porepressure
value<<<{"description": "Value of the BC"}>>> = 1E6
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
[]
[constant_outer_porepressure]
type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = porepressure
value<<<{"description": "Value of the BC"}>>> = 0
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = rmax
[]
[injected_tracer]
type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = tracer_concentration
value<<<{"description": "Value of the BC"}>>> = 0.5
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
[]
[]
[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
[the_simple_fluid]
type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../source/fluidproperties/SimpleFluidProperties.html"}>>>
bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2E9
viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1.0E-3
density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000.0
[]
[]
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[porosity]
type = PorousFlowPorosity<<<{"description": "This Material calculates the porosity PorousFlow simulations", "href": "../../source/materials/PorousFlowPorosity.html"}>>>
porosity_zero<<<{"description": "The porosity at zero volumetric strain and reference temperature and reference effective porepressure and reference chemistry. This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)"}>>> = 0.1
[]
[permeability_aquifer]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquifer
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[permeability_caps]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = caps
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-15 0 0 0 1E-15 0 0 0 1E-16'
[]
[]
[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
active<<<{"description": "If specified only the blocks named will be visited and made active"}>>> = basic
[basic]
type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = ' asm lu NONZERO 2'
[]
[preferred_but_might_not_be_installed]
type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = ' lu mumps'
[]
[]
[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_rel_tol = 1E-14
[]
[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[]The only other changes are that the PorousFlowPorosity Material is needed, not the PorousFlowConstantBiotModulus Material:
[Materials]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
If multiply_by_density = true then the arguments in Page 02 concerning the magnitude of the residual, and hence the size of the nonlinear absolute-tolerance must be modified. Specifically, the fluid residuals now get multiplied by the density (which is the default and most common use-case in PorousFlow). So in this problem, the fluid residual is approximately (since the fluid density is approximately 1000kg.m). The nl_abs_tol is therefore
[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_rel_tol = 1E-14
[]An animation of the results is shown in Figure 1.
Figure 1: Darcy flow vectors and tracer evolution in the borehole-aquifer-caprock system.
As you gain experience with PorousFlow, you will realise that this simulation contains several defects.
Boundary conditions
The most important (potential) problem lies with the boundary conditions. Physically they are saying "add or remove fluid to keep the porepressure or tracer concentration fixed". In more complex simulations this can be completely disastrous. There may not be the correct fluid component and fluid phase present to remove, and this causes extremely poor nonlinear convergence. It is really important to understand this fully.
PorousFlowFullySaturated associates the variable with the mass-balance equation for , and the variable with the mass-balance equation for the final mass fraction (). (This is not obvious unless you read the documentation or the code. In simulations that don't use Actions it will be obvious because you will have to write the input-file Kernels block yourself.) Hence the BCs that have variable = porepressure:
[constant_injection_porepressure]
type = DirichletBC
variable = porepressure
value = 1E6
boundary = injection_area
[]
[constant_outer_porepressure]
type = DirichletBC
variable = porepressure
value = 0
boundary = rmax
[]are saying "add or remove the final mass fraction () to keep the porepressure fixed" (physically this mass fraction is probably the water mass fraction); while the BCs that have variable = tracer_concentration:
[injected_tracer]
type = DirichletBC
variable = tracer_concentration
value = 0.5
boundary = injection_area
are saying "add or remove the tracer so that its concentration remains fixed". In complex multi-phase, multi-component models this can be a real "gotcha" and a PorousFlowSink is recommended.
Lack of numerical stabilization
Broadly speaking, the animation shows the expected behaviour, but what are those "stripes" of low and high concentration towards the end of the simulation? This is caused by the lack of numerical stabilization. As mentioned in the introduction to this Page, stabilization is standard in most PorousFlow simulations, but is not used in the above simulation. In fact, if you run the simulation yourself, you'll see that tracer_concentration does not even lie within its expected range of !
The "stripes" are removed by numerical stabilization, but this comes at a cost. PorousFlow offers two types of numerical stabilization. Full-upwinding is excellent at eliminating overshoots and undershoots, and the convergence speed is also excellent, but it introduces extra diffusion, so that sharp fronts are not maintained! On the other hand, KT stabilization also eliminates overshoots and undershoots and introduces very little or zero extra diffusion so that sharp fronts are maintained as well as without any numerical stabilization, but its convergence speed and runtime are always greater. Further details regarding the diffusion of fronts may be found in the numerical diffusion page.
Numerical stabilization can be easily included using the PorousFlowFullySaturated Action. For instance, KT stabilization may be easily introduced:
[PorousFlowFullySaturated]
porepressure = porepressure
coupling_type = Hydro
gravity = '0 0 0'
fp = the_simple_fluid
mass_fraction_vars = tracer_concentration
stabilization = KT
flux_limiter_type = superbee
Users are strongly encouraged to experiment with KT stabilization and report their findings back to the moose-users google group, since this stabilization is still at the development stage.