Porosity based on STATE (MaterialProperty)- cant figure out, what did i do wrong

[Traiwit]

Hi PorousFlow experts,

I need help with my Porosity class (where the porosity is updated based on state value) where STATE is another MaterialProperty

let's look at my simple input file below

## mpiexec -n 4 ~/moose_beh/moose_beh-opt -i groundwater_test.i --n-threads=4 --t

gravity = '0 0 -9.81'
[Mesh]
  [basic_mesh]
    type = GeneratedMeshGenerator
    dim = 3
    xmin = -200
    xmax = 200
    nx = 10
    ymin = -50
    ymax = 50
    ny = 5
    zmin = -200
    zmax = -100
    nz = 5
  []

  [block_2]
    type = SubdomainBoundingBoxGenerator
    input = basic_mesh
    block_id = 2
    block_name = lower_aquifer
    bottom_left = '150 -50 -200'
    top_right = '200 50 -100'
  []
[]

[Problem]
  type = 'FEProblem'   # (see https://mooseframework.inl.gov/source/problems/FEProblem.html)
  
  ### checking parameters:
  boundary_restricted_elem_integrity_check = true   # default: true
  boundary_restricted_node_integrity_check = true   # default: true
  check_uo_aux_state = false                        # default: false
  error_on_jacobian_nonzero_reallocation = false    # default: false
  fv_bcs_integrity_check = true                     # default: true
  kernel_coverage_check = true                      # default: true
  material_coverage_check = true                    # default: true
  material_dependency_check = false                  # default: true
  skip_nl_system_check = false                      # default: false
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
    [dictator]
      type = PorousFlowDictator
      porous_flow_vars = pp
      number_fluid_phases = 1
      number_fluid_components = 1
      execute_on = 'TIMESTEP_END'   # the default
    []
  []

[Variables]
  [pp]
  []
[]

[ICs]
  [pp]
    type = FunctionIC
    variable = pp
    function = fcn_pwp_initialCondition
  []
[]

[BCs]
  [pp]
    type = FunctionDirichletBC
    variable = pp
    function = fcn_pwp_initialCondition
    boundary = 'left'
  []

  # [bc_drain]
  #   type = PorousFlowBEH_Drainage
  #   boundary = 'right'
  #   variable = pp
  #   # char_elem_length_n = charElLength_M0
  #   fluid_phase = 0
  #   pt_vals =     '0 1e9'
  #   multipliers = '0 1e9'
  #   PT_shift = 0.
  #   flux_function = 1e-1
  #   use_mobility = true   
  #   use_relperm = false   
  # []
[]

[Functions]
  [fcn_pwp_initialCondition]
    type = ParsedFunction
    expression = '9810.*(-100-z)'
  []
[]

[AuxVariables]
  [porosity_M0]
    type = MooseVariableConstMonomial
  []
[]

[AuxKernels]
  [porosity_M0]
    type = PorousFlowPropertyAux
    variable = porosity_M0
    property = porosity
    execute_on = 'INITIAL TIMESTEP_BEGIN TIMESTEP_END'
  []
[]



[Kernels]
    [fluid_mass_dot]
      type = PorousFlowMassTimeDerivative
      variable = pp
      fluid_component = 0
      multiply_by_density = true
    []
    
    [fluid_flux_fullupwind]
      type = PorousFlowAdvectiveFlux
      variable = pp
      gravity = '${gravity}'
      fluid_component = 0
    []
  []

[FluidProperties]
  [simple_fluid]
    type = SimpleFluidProperties
    thermal_expansion = 0
    bulk_modulus = 1e9
    density0 = 1000
  []
[]

[UserObjects]
    [uobj_saturation]
      type = PorousFlowCapillaryPressureVG
      alpha = 1e-6
      m = 0.6
      log_extension = True
      sat_lr = 0.1
    []
[]

[Materials]
  [mat_groundwater]
      type = PorousFlowSingleComponentFluid
      fp = simple_fluid
      phase = 0
      pressure_unit = Pa          # options: Pa (default), MPa
      temperature_unit = Kelvin   # options: Kelvin (default), Celsius
      time_unit = years           # options: seconds (default), hours, days, years
      compute_density_and_viscosity = true
  []
  [mat_temperature]
    type = PorousFlowTemperature
  []
  [mat_massfrac]
    type = PorousFlowMassFraction
  []

  [mat_saturation]
    type = PorousFlow1PhaseP
    porepressure = pp
    capillary_pressure = uobj_saturation
  []
  
  [porosity_block0]
    type = PorousFlowPorosityConst
    porosity = 0.07
    block = 0
  []
  [porosity_block2_1]
    type = PorousFlowPorosityConst
    porosity = 0.075
    block = 2
  []

  # [PorousFlowPorosityFromCSV]
  #   type = PorousFlowPorosityFromCSV
  #   state = state
  #   state_values =      '0    1       2    3'
  #   poro_state_values = '0.07 0.07  0.075 0.1'
  # []

    ## specify state here
  [state1]
    type = GenericFunctionMaterial
    block = 0
    prop_names = state
    prop_values = 1
  [../]
  [state2]
    type = GenericFunctionMaterial
    block = 2
    prop_names = state
    prop_values = 2
  [../]

  [biot_mod]
    type = PorousFlowConstantBiotModulus
    fluid_bulk_modulus = 2E9
    biot_coefficient = 1.0
  []
  [permeability_aquifers]
    type = PorousFlowPermeabilityConst
    permeability = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
  []
  [mat_relative_permeability]
      type = PorousFlowRelativePermeabilityConst
      phase = 0
      kr = 1.0
  []
[]

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

[Executioner]
  type = Transient
  solve_type = Newton
  use_pre_SMO_residual = false
  dt = 0.0005
  end_time = 0.01 
  nl_abs_tol = 1E-5
  nl_rel_tol = 1E-20
[]

[Outputs]
  [exodus]
    type = Exodus
    file_base = test_poro_change
    output_material_properties = true
  []
[]

so as you can see, there are 2 ways of assigning the porosity:

1. MOOSE standard way

  [porosity_block0]
    type = PorousFlowPorosityConst
    porosity = 0.07
    block = 0
  []
  [porosity_block2_1]
    type = PorousFlowPorosityConst
    porosity = 0.075
    block = 2
  []

2. from my class

 # [PorousFlowPorosityFromCSV]
  #   type = PorousFlowPorosityFromCSV
  #   state = state
  #   state_values =      '0    1       2    3'
  #   poro_state_values = '0.07 0.07  0.075 0.1'
  # []

    ## specify state here
  [state1]
    type = GenericFunctionMaterial
    block = 0
    prop_names = state
    prop_values = 1
  [../]
  [state2]
    type = GenericFunctionMaterial
    block = 2
    prop_names = state
    prop_values = 2
  [../]

the porosity showing on Paraview are exactly the same, but the pp are drastically different (only when PorousFlowMassTimeDerivative is activated)

pp with assign porosity from MOOSE default method

PP when assign porosity from my class

and below is my class .C

i couldnt figure out, what did i do wrong? also you can see the density are a lil bit different too

//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "PorousFlowPorosityFromCSV.h"

registerMooseObject("moose_behApp", PorousFlowPorosityFromCSV);

InputParameters
PorousFlowPorosityFromCSV::validParams()
{
  InputParameters params = PorousFlowPorosityBase::validParams();
  params.addParam<UserObjectName>("read_prop_user_object",
                                  "The ElementReadPropertyFile GeneralUserObject to read element-specific property values from file");
  params.addClassDescription("This Material calculates the porosity assuming it is constant");
  params.addParam<MaterialPropertyName>("state", 0, "state");
  params.addParam<std::vector<int>>("state_values", std::vector<int>(), "List of state values");
  params.addParam<std::vector<Real>>("poro_state_values", std::vector<Real>(), "List of porosity values corresponding to the state values");
  return params;
}

PorousFlowPorosityFromCSV::PorousFlowPorosityFromCSV(const InputParameters & parameters)
  : PorousFlowPorosityBase(parameters),
    _read_prop_user_object(isParamValid("read_prop_user_object")
                               ? &getUserObject<ElementPropertyReadFile>("read_prop_user_object")
                               : nullptr),
    _state(getMaterialProperty<Real>("state")),
    _stateValues(getParam<std::vector<int>>("state_values")),
    _poroStateValues(getParam<std::vector<Real>>("poro_state_values"))
{
    // Ensure the state and poro_state vectors are of the same size
    if (_stateValues.size() != _poroStateValues.size()) {
        mooseError("state_values and poro_state_values must have the same length");
    }

    // Populate the state to porosity map
    for (size_t i = 0; i < _stateValues.size(); ++i) {
        _stateToPoroMap[_stateValues[i]] = _poroStateValues[i];
    }
}

void
PorousFlowPorosityFromCSV::initQpStatefulProperties()
{
  // Calculate Porosity
  if (_stateToPoroMap.find(_state[_qp]) != _stateToPoroMap.end()) {
    // If the state exists in the map, use the corresponding porosity value
    Real poroValue = _stateToPoroMap[_state[_qp]];
    _porosity[_qp] = poroValue;
  } else {
    if (_read_prop_user_object) {
      _porosity[_qp] = _read_prop_user_object->getData(_current_elem, 0);
    } else {
      mooseWarning("No user object specified for reading property values, and state value not found in the map. Porosity will remain undefined.");
    }
  }
}

void
PorousFlowPorosityFromCSV::computeQpProperties()
{
  initQpStatefulProperties();

    // The derivatives are zero for all time
    (*_dporosity_dvar)[_qp].assign(_num_var, 0.0);
    (*_dporosity_dgradvar)[_qp].assign(_num_var, RealGradient());
}

thank you!

Hi @cpgr I tag you since you are the expert here 🥇

quick note:

it seems that the problem only occurs when i assign porosity based on state (like below)

state_values = '0 1 2 3'
poro_state_values = '0.07 0.07 0.075 0.1'

varied spatial porosity from .csv file does not give this weird pp/results

[Traiwit]

interesting finding

so i let it run a bit, the pp seems to go back to where it should be

not even sure why, since

state_values = '0 1 2 3'
poro_state_values = '0.07 0.07 0.075 0.1'

are assigned since INITIAL similar to MOOSE default way

[Traiwit]

I have been thinking of your problem and believe it would be better to make an AuxVariable porosity that gets set based on a state AuxVariable. You can probably achieve what you want using a FunctionAux with a PiecewiseConstant function? The advantage of this is that the porosity AuxVariable can be used with the existing PorousFlowPorosityConst material to provide a spatially heterogeneous porosity based on the state variable.

thank you for your reply @cpgr

if I do that, in PorousFlowMassTimeDerivative do i just grab an updated porosity from the AuxVariable instead of MaterialProperty?

[cpgr]

No - the good part of using an AuxVar is that you don't have to anything other than

[porosity]
  type = PorousFlowPorosityConst
  porosity = porosity_var_name
[]

That material will then grab the latest value of the AuxVar and all of the old values should get propagated correctly!

[Traiwit]

@cpgr awesome, thank you!

i managed to update porosity_nodal using moose default method, but i dont really like the results

the way i did it: i added another set of PorousFlowPorosityConst with at_nodes = true, i checked porosity_nodal now is updated correctly (also double checked with Paraview) - please see below

what i dont like: changing porosity between time_step seems to drain out the water (gives sudden negative pp) and if i let it run, it will reflood the void and go back to normal with time (idk kinda make sense right? when the pore space suddenly get much bigger but amount of fluid remain the same)

first time step:

at some point:

Now i have 2 choice for my larger scale hydro-mining simulation: use constant porosity - ignore the change porosity in excavated zone

2: disable PorousFlowMassTimeDerivative at the excavated block

note excavated zone: porosity = 1 state = 1

    [fluid_mass_dot]
      type = PorousFlowMassTimeDerivative
      variable = pp
      fluid_component = 0
      multiply_by_density = true
      block = 0
    []

 [porosity_block0]
    type = PorousFlowPorosityConst
    porosity = 0.07
    block = 0
  []
  [porosity_block2_1]
    type = PorousFlowPorosityConst
    porosity = 0.071
    block = 2
  []
  [porosity_block2_2]
    type = PorousFlowPorosityConst
    porosity = 0.85
    block = 2
    enable = false
  []
  [porosity_block2_3]
    type = PorousFlowPorosityConst
    porosity = 1
    block = 2
    enable = false
  []

  [porosity_block0_node]
    type = PorousFlowPorosityConst
    porosity = 0.07
    block = 0
    at_nodes = true
  []
  [porosity_block2_1_node]
    type = PorousFlowPorosityConst
    porosity = 0.071
    block = 2
    at_nodes = true
  []
  [porosity_block2_2_node]
    type = PorousFlowPorosityConst
    porosity = 0.85
    block = 2
    enable = false
    at_nodes = true
  []
  [porosity_block2_3_node]
    type = PorousFlowPorosityConst
    porosity = 1
    block = 2
    enable = false
    at_nodes = true
  []

  [control_1]
    type = ConditionalFunctionEnableControl
    conditional_function = 'if(t > 0.0005 & t<=0.001 ,1,0)'
    enable_objects = 'Materials::porosity_block2_2 Materials::porosity_block2_2_node'
    disable_objects = 'Materials::porosity_block2_1 Materials::porosity_block2_1_node'
    execute_on = 'INITIAL TIMESTEP_BEGIN TIMESTEP_END'
  [../]
  [control_2]
    type = ConditionalFunctionEnableControl
    conditional_function = 'if(t>0.001,1,0)'
    enable_objects = 'Materials::porosity_block2_3 Materials::porosity_block2_3_node'
    disable_objects = 'Materials::porosity_block2_2 Materials::porosity_block2_2_node'
    execute_on = 'INITIAL TIMESTEP_BEGIN TIMESTEP_END'
  [../]

[cpgr]

what i dont like: changing porosity between time_step seems to drain out the water (gives sudden negative pp) and if i let it run, it will reflood the void and go back to normal with time (idk kinda make sense right? when the pore space suddenly get much bigger but amount of fluid remain the same)

I think this is probably expected - making the porosity larger can mean that the pressure drops and could lead to desaturation (negative pressure in this case by convention). You can force the rock to remain saturated using the fully saturated version of the pressure-saturation material, but the pressure will still drop as porosity increases so that the fluid density decreases to keep the rock saturated. I guess it depends what behaviour you are trying to model?