Deactivate/delete elements based on material property

[Traiwit]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hi guys,

just wondering is there a way to deactivate/delete element based on material property?

for example, if porosity of this element = 1, then don't solve/ignore this element in the system (mimic the excavation)

I checked CoupledVarThresholdElementSubdomainModifier this is more like re-grouping blocks to me, my simulation is element-based, all material properties are imported via .csv, I kinda don't want to deal with block (element-set) system.

Thank you!
Traiwit

[GiudGiud]

Hello

No it hasnt been added.
You can add a new object, or you can convert a material property to an auxiliary variable with an MaterialRealAux for example then use the current way
CoupledVarThresholdElementSubdomainModifier is exactly how we do this currently. You re-assign the elements you do not care to solve on to a block which does not have any kernels / equations defined on

Guillaume

[GiudGiud]

Hello

A few questions,

• you are removing elements from the calculation right?
  Should you be moving a boundary at the same time to keep imposing a boundary condition where you are removing?

• What does the variable block look like? Did you make sure to block restrict to block 1 and 2?

• Why execute_on = 'TIMESTEP_BEGIN TIMESTEP_END' ? Do you have something executing in between (use Debug/show_execution_order to check maybe)? Because if not you are executing twice in a row (end of a time step, then beginning of the next one).

• The nonlinear residual even without the UO is not going down that low. Is this as far as you can get it to converge?

• Even without the Uo the linear solve is reporting convergence difficulties, have you tried using a different preconditoner?

[Traiwit]

Hi @GiudGiud

you are removing elements from the calculation right?
Should you be moving a boundary at the same time to keep imposing a boundary condition where you are removing?

The boundaries are predefined via Abaqus mesh, which exactly wrapping around those elements at each step:

What does the variable block look like? Did you make sure to block restrict to block 1 and 2?

I wasn't aware of this, I will try this and report back:

[Variables]
  [porepressure_L1]
    type = MooseVariable
    order = FIRST
    family = LAGRANGE
    block = '1 2'
  []
[]

Why execute_on = 'TIMESTEP_BEGIN TIMESTEP_END' ? Do you have something executing in between (use Debug/show_execution_order to check maybe)? Because if not you are executing twice in a row (end of a time step, then beginning of the next one).

TIMESTEP_BEGIN should be enough, since STATE will be updated at TIMESTEP_BEGIN as well (read from csv - element based)

The nonlinear residual even without the UO is not going down that low. Is this as far as you can get it to converge?
Even without the Uo the linear solve is reporting convergence difficulties, have you tried using a different preconditoner?

yep first step will mostly be taken care by nl_rel_tol, but later it will go as low as 1e4 to 1e5, this is as low as i can push it down, the system is very nonlinear due to (partially saturation and discontinuity). Below is my solver setup

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient   # always use 'Transient' (not 'Steady'), see [Kernels]
  solve_type = NEWTON
  petsc_options = '-snes_converged_reason'
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre    boomeramg'
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e5
  nl_max_its = 20
  nl_forced_its = 3
  l_tol = 1e-6
  l_max_its = 200
  # line_search = 'none'
  ###
  ### time stepping:
  # dt = 1.0
  dtmin = 0.083
  dtmax = 1.0
  start_time = 1.75
  end_time = 3
  [TimeStepper]
    ### option-2: define time steps here, directly
    type = TimeSequenceStepper
    time_sequence = '1.000 1.500 1.750 1.833 2.000 2.083 2.167 2.250 2.333 2.417 2.500 . . . . . very long
  []
[]

I can also send you the full input file if you'd like, but it's quite long like 5000+ lines

Thank you for you help!

Kind regards,
Traiwit

[GiudGiud]

The boundaries are predefined via Abaqus mesh, which exactly wrapping around those elements at each step:

I m not sure I understand this properly as it seems unlikely that Abaqus knows which elements are getting removed (moved to an inactive block) beforehand, and the boundary at the edge of the removal could not be defined this way?
What is the boundary condition at the boundary near removed (moved to an inactive block) elements?

[Traiwit]

Hi @GiudGiud, so we predefine the excavation and drain surface for each timestep in Abaqus mesh. The image below shows the surfaces not the blocks

then we use multiple BCs + Controls to change the drain surface at each timestep, for example:

  [bc_drain0005_UG]
    type = PorousFlowPiecewiseLinearSinkBEH
    variable = porepressure_L1
    boundary = 'DRAINSRFUG_F005_Y2019_M12'
    pt_vals = '0 1e9'
    multipliers = '0 1e9'
    flux_function = 1
    #charElLength = charElLength_M0
    PT_shift = 0.
    fluid_phase = 0
    use_mobility = true    # default is false
    use_relperm = false     # default is false
  []
  
  ### --- TOPO-RECHARGE AND DRAINAGE FOR 0006: F006_Y2020_M01
  ### NOT DEFINING [bc_topol0006_LGFlux]
  ### NOT DEFINING [bc_topol0006_Rainfall]
  ### NOT DEFINING [bc_topol0006_RiverLake]
  ### NOT DEFINING [bc_drain0006_Pit]
  [bc_drain0006_UG]
    type = PorousFlowPiecewiseLinearSinkBEH
    variable = porepressure_L1
    boundary = 'DRAINSRFUG_F006_Y2020_M01'
    pt_vals = '0 1e9'
    multipliers = '0 1e9'
    flux_function = 1
    #charElLength = charElLength_M0
    PT_shift = 0.
    fluid_phase = 0
    use_mobility = true    # default is false
    use_relperm = false     # default is false
  []
  
  ### --- TOPO-RECHARGE AND DRAINAGE FOR 0007: F007_Y2020_M02
  ### NOT DEFINING [bc_topol0007_LGFlux]
  ### NOT DEFINING [bc_topol0007_Rainfall]
  ### NOT DEFINING [bc_topol0007_RiverLake]
  ### NOT DEFINING [bc_drain0007_Pit]
  [bc_drain0007_UG]
    type = PorousFlowPiecewiseLinearSinkBEH
    variable = porepressure_L1
    boundary = 'DRAINSRFUG_F007_Y2020_M02'
    pt_vals = '0 1e9'
    multipliers = '0 1e9'
    flux_function = 1
    #charElLength = charElLength_M0
    PT_shift = 0.
    fluid_phase = 0
    use_mobility = true    # default is false
    use_relperm = false     # default is false
  []

  [controls0005]
    type = ConditionalFunctionEnableControl
    enable_objects = '    BCs::bc_drain0005_UG'
    # conditional_function = 'if(t>=2.000 & t<2.083,1,0)'
    # execute_on = 'INITIAL TIMESTEP_BEGIN TIMESTEP_END'
    conditional_function = 'if(t>1.833 & t<=2.000,1,0)'
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
  [controls0006]
    type = ConditionalFunctionEnableControl
    enable_objects = '    BCs::bc_drain0006_UG'
    # conditional_function = 'if(t>=2.083 & t<2.167,1,0)'
    # execute_on = 'INITIAL TIMESTEP_BEGIN TIMESTEP_END'
    conditional_function = 'if(t>2.000 & t<=2.083,1,0)'
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
  [controls0007]
    type = ConditionalFunctionEnableControl
    enable_objects = '    BCs::bc_drain0007_UG'
    # conditional_function = 'if(t>=2.167 & t<2.250,1,0)'
    # execute_on = 'INITIAL TIMESTEP_BEGIN TIMESTEP_END'
    conditional_function = 'if(t>2.083 & t<=2.167,1,0)'
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []

The BC is exactly like PorousFlowPiecewiseLinearSink when charElLength = charElLength_M0 isnt provided.

we want to do inactive elements because we think it might help with convergence + reduce the RAM requirement a bit.

The simulation is working fine without inactive elements but we are just not sure what material property to assign to those elements to mimic the inactive behaviour.

[GiudGiud]

The simulation is working fine without inactive elements but we are just not sure what material property to assign to those elements to mimic the inactive behaviour.

@cpgr do you know