PorousFlowPolyLineSink not working with plasticity

hsheldon commented 5 years ago

Rationale

PorousFlowPolyLineSink should work with plasticity @WilkAndy

Description

The following input file seg faults with git commit 7bb99349a8, but was working with git commit 7eb99a048b

# Example: Injection into a uniform aquifer 10 x 10 x 5 km
# Drucker-Prager deformation
# Darcy flow
# Using PorousFlow module
# Small example to test seg fault

[Mesh]
  type = GeneratedMesh
  dim = 3
  xmin = 0
  xmax = 1e4
  ymin = 0
  ymax = 1e4
  zmin = -5e3
  zmax = 0
  nx = 6
  ny = 6
  nz = 6
[]

[GlobalParams]
  block = '0'

  PorousFlowDictator = dictator
  gravity = '0 0 -9.81'
  biot_coefficient = 1.0
  solid_bulk = 1.0 # Required but irrelevant when biot_coefficient is unity

  displacements = 'disp_x disp_y disp_z'

  strain_at_nearest_qp = true

  # used by TensorMechanicsPlasticDruckerPragerHyperbolic
  yield_function_tolerance = 1
  internal_constraint_tolerance = 1E-6
  mc_interpolation_scheme = inner_tip # outer_tip, inner_tip or inner_edge
  use_custom_returnMap = true

  # Used in DruckerPrager materials:
  max_NR_iterations = 1000
  smoothing_tol = 1E5
  yield_function_tol = 1E-3
  perform_finite_strain_rotations = false
[]

[Modules]
  [./FluidProperties]
    [./simple_fluid]
      type = SimpleFluidProperties
      thermal_expansion = 0 # Set this to zero if you are not doing a thermal simualtion
      bulk_modulus = 2E9
      density0 = 1000
      viscosity = 5E-4
    [../]
  [../]
[]

[PorousFlowFullySaturated]
  coupling_type = HydroMechanical
  porepressure = pp
  dictator_name = dictator
  fp = simple_fluid
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
  [./pp]
    scaling = 1E6
    [./InitialCondition]
      type = FunctionIC
      function = ini_pp
    [../]
  [../]
[]

[Functions]
  [./ini_stress]
    type = ParsedFunction
    vars = 'g rho_s rho_f poro' # gravity, solid density, fluid density, porosity
    vals = '-9.81 2350 1000 0.1'
    value = '-g*z*(rho_s-rho_f)*(1.0-poro)' # initial effective stress that should result from weight force
  [../]

  # Initial fluid pressure
  [./ini_pp]
    type = ParsedFunction
    vars = 'g rho_s rho_f'
    vals = '-9.81 2350 1000' # gravity, solid density, fluid density
    value = 'g*z*rho_f+1E5' 
  [../]
[]

[BCs]
  [./p_top]
    type = FunctionPresetBC
    variable = pp
    boundary = front
    function = ini_pp
  [../]

  [./xsides]
    type = PresetBC
    variable = disp_x
    boundary = 'left right'
    value = '0'
  [../]
  [./ysides]
    type = PresetBC
    variable = disp_y
    boundary = 'top bottom'
    value = '0'
  [../]
  [./z_sides]
    type = PresetBC
    variable = disp_z
    boundary = 'back front'
    value = '0'
  [../]
[]

[AuxVariables]
[]

[AuxKernels]
[]

[Postprocessors]
[]

[UserObjects]
  [./pls_total_outflow_mass]
    type = PorousFlowSumQuantity
  [../]

  # Cohesion
  [./mc_coh]
    type = TensorMechanicsHardeningConstant
    value = 6.0E6
  [../]

  # Friction angle
  [./mc_phi]
    type = TensorMechanicsHardeningConstant
    value = 35.0
    convert_to_radians = true
  [../]

  # Dilation angle
  [./mc_psi]
    type = TensorMechanicsHardeningConstant
    value = 2
    convert_to_radians = true
  [../]

  # Drucker-Prager objects
  [./dp]
    type = TensorMechanicsPlasticDruckerPragerHyperbolic
    mc_cohesion = mc_coh
    mc_friction_angle = mc_phi
    mc_dilation_angle = mc_psi
    smoother = 1E6
  [../]

  # Tensile strength
  [./tens]
    type = TensorMechanicsHardeningConstant
    value = 3.0E6
  [../]

  # Compressive strength (cap on yield envelope)
  [./compr_all]
    type = TensorMechanicsHardeningConstant
    value = 1E10
  [../]
[]

[Materials]  
  [./strain]
    type = ComputeIncrementalSmallStrain
    eigenstrain_names = eigenstrain_all
  [../]

  [./eigenstrain_all]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = 'ini_stress 0 0  0 ini_stress 0  0 0 ini_stress'
    eigenstrain_name = eigenstrain_all
  [../]

  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    bulk_modulus = 3.3333E9
    shear_modulus = 2.5E9
  [../]

  [./dp_mat]
    type = CappedDruckerPragerStressUpdate
    DP_model = 'dp'
    tensile_strength = 'tens'
    compressive_strength = 'compr_all'
    tip_smoother = 3.0E6
  [../]

  [./stress]
    type = ComputeMultipleInelasticStress
    inelastic_models = dp_mat
  [../]

  # Permeability
  [./permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-13 0 0  0 1E-13 0  0 0 1E-13'
  [../]

  # Porosity 
  [./porosity_nodal]
    type = PorousFlowPorosity
    at_nodes = true
    porosity_zero = 0.1
    mechanical = true
    fluid = true
  [../]

  [./porosity_qp]
    type = PorousFlowPorosity
    porosity_zero = 0.1
    mechanical = true
    fluid = true
  [../]

 # Density of saturated rock
  [./density]
    type = PorousFlowTotalGravitationalDensityFullySaturatedFromPorosity
    rho_s = 2350
  [../]
[]

[DiracKernels]
  [./pls]
    type = PorousFlowPolyLineSink
    variable = pp
    SumQuantityUO = pls_total_outflow_mass
    point_file = two_nodes2.bh
    function_of = pressure
    fluid_phase = 0   
    p_or_t_vals = '0 1E7'
    fluxes = '-1.59 -1.59'
  [../]
[]

[Preconditioning]
  [./andy]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  solve_type = Newton
  type = Transient
  line_search = bt

  dt = 2e5
  end_time = 1e6

  l_tol = 1E-5
  nl_abs_tol = 1E0
  nl_rel_tol = 1E-8
  l_max_its = 200
  nl_max_its = 30

  petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
  petsc_options_value = ' asm      2              lu            gmres     200'
[]

[Outputs]
  file_base = uniform_genmesh_inject_SEGFAULT
  exodus = true
  csv = true
  execute_on = 'initial timestep_end'
[]

Impact

Necessary for users wishing to couple plasticity with PorousFlowPolyLineSink

hsheldon commented 5 years ago

Here's the .bh file for the above example:

1 5000.0 5000.0 -1400.0
1 5000.0 5000.0 -1500.0