Chessboard pattern in pore pressure

[jmeier]

Dear Moose Community,

To get a little more insight what is going on in #30828, I did some tinkering with one of the porous flow tests: modules/porous_flow/test/tests/sinks/s06.i. Up to my understanding this is a pure porous-flow test with no solid mechanics involved.

My context is: I have to build a hydro-mechanical coupled Moose model using TET10 elements.

Originally, this test consists of one HEX8 element (nx=ny=nz=1). I modified to have nx=ny=nz=5 and second order elements (HEX27 or TET10). Furthermore, I modified some of the material constants (to get more physical values), scaled dt (and end_time) and reduced tolerances a bit. But nothing what should change the overall behavior. Please find the resulting input file below.

Switching between HEX27 and TET10 leads to some differences starting from time step 1. Please see the image below. The result for the HEX27 seem to be plausible. But it looks as if for TET10 the porous flow sink is not applied to the nodes at the corners of the elements but to the nodes at the midpoints:

1. Is it correct that the behavior for TET10 is not plausible?
2. Is there something special to be done to improve the system response for TET10?

input file

# apply a half-cubic sink flux and observe the correct behavior
[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 5
  ny = 5
  nz = 5
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
  zmin = 0
  zmax = 2
  elem_type = HEX27 #TET10
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 0.0000000005 #1.1
  []
[]

[Variables]
  [pp]
    family = LAGRANGE
    order = SECOND
  []
[]

[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    function = 'x * (y + 1)'
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 1.3e9
    density0 = 1100
    thermal_expansion = 0
    viscosity = 1.1e-3
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-5 0 0 0 1E-5 0 0 0 1E-5'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
[]

[AuxVariables]
  [flux_out]
    family = LAGRANGE
    order = SECOND
  []
[]

[BCs]
  [flux]
    type = PorousFlowHalfCubicSink
    boundary = 'left right'
    max = 2
    cutoff = -0.8
    center = 0.9
    variable = pp
    use_mobility = false
    use_relperm = false
    fluid_phase = 0
    flux_function = 3
    save_in = flux_out
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -snes_max_it -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' gmres     asm      lu           10000        NONZERO                   2'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = ' lu       mumps'


  dt = 2000
  end_time = 60000

  nl_rel_tol = 1E-9 #12
  nl_abs_tol = 1E-9 #12
[]

[Outputs]
  exodus = true
  [console]
    type = Console
    execute_on = 'nonlinear linear'
    time_step_interval = 5
  []
[]

[Debug]
  # check_jacobian = true
  show_top_residuals = 0
  show_var_residual_norms = true
[]

[GiudGiud]

Hello

Given that hex20 looks weird too and tet14 segfaults (oob in the material properties), I think this is a problem with porous flow and higher order variables. Maybe the nodal mass lumping

@jmeier can we simplify this further? Maybe using less material properties, using standard BCs and kernels.
What is the equivalent of a PorousFlowHalfCubicSink? It's just a flux right (sink? so negative)

I remember @lindsayad was interested in tet10 oscillations with constraints so here's a ping.

[jmeier]

The model given in #30698 also uses TET10 and shows plausible results in case Quadrature=SIMPSON is not used.

[cpgr]

Hmmm, I've been looking at this, and am wondering if these affects are due to the really simplified physics included. Some of these tests are really contrived to test one thing. In the above examples, there is no fluid flux, only a mass balance. And the mass balance is lumped to the nodes, which perhaps leads to the checkerboarding.

I've played with a more realistic case of flux into a model (attached below) and everything seems to work with first or second order variables on hexes, tet10's and hex27' with an integrated BC.

and with the meshes (hex top, tet10 middle, hex27 bottom)

[Mesh]
  type = GeneratedMesh
  dim = 3
  nx = 50
  ny = 3
  nz = 3
  xmin = 0
  xmax = 100
  ymax = 10
  zmax = 10
  # elem_type = TET10
  elem_type = HEX27
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pp]
    initial_condition = 2E6
    family = LAGRANGE
    order = SECOND
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pp
  []
  [flux]
    type = PorousFlowAdvectiveFlux
    variable = pp
    gravity = '0 0 0'
    fluid_component = 0
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e-7
  []
[]

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2e9
    density0 = 1000
    thermal_expansion = 0
    viscosity = 1e-3
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
  []
  [ppss]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [simple_fluid]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid
    phase = 0
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
  []
  [relperm]
    type = PorousFlowRelativePermeabilityCorey
    n = 0
    phase = 0
  []
[]

[BCs]
  # [left]
  #   type = DirichletBC
  #   boundary = left
  #   value = 3E6
  #   variable = pp
  # []
  [left]
    type = PorousFlowSink
    boundary = left
    flux_function = -1e-4
    variable = pp
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type'
    petsc_options_value = 'bcgs bjacobi'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E3
  end_time = 1E4
[]

[Postprocessors]
  [p000]
    type = PointValue
    variable = pp
    point = '0 0 0'
    execute_on = 'initial timestep_end'
  []
  [p010]
    type = PointValue
    variable = pp
    point = '10 0 0'
    execute_on = 'initial timestep_end'
  []
  [p020]
    type = PointValue
    variable = pp
    point = '20 0 0'
    execute_on = 'initial timestep_end'
  []
  [p030]
    type = PointValue
    variable = pp
    point = '30 0 0'
    execute_on = 'initial timestep_end'
  []
  [p040]
    type = PointValue
    variable = pp
    point = '40 0 0'
    execute_on = 'initial timestep_end'
  []
  [p050]
    type = PointValue
    variable = pp
    point = '50 0 0'
    execute_on = 'initial timestep_end'
  []
  [p060]
    type = PointValue
    variable = pp
    point = '60 0 0'
    execute_on = 'initial timestep_end'
  []
  [p070]
    type = PointValue
    variable = pp
    point = '70 0 0'
    execute_on = 'initial timestep_end'
  []
  [p080]
    type = PointValue
    variable = pp
    point = '80 0 0'
    execute_on = 'initial timestep_end'
  []
  [p090]
    type = PointValue
    variable = pp
    point = '90 0 0'
    execute_on = 'initial timestep_end'
  []
  [p100]
    type = PointValue
    variable = pp
    point = '100 0 0'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  print_linear_residuals = false
  exodus = true
[]

[cpgr]

So perhaps the problem is something else in your more realistic model. I'd like to help out there if possible.

[jmeier]

As a side-note: The model given in #30828 is also more realistic and also shows numerical instabilities.

[cpgr]

I'll take a look at this tomorrow my time