search

idaholab / moose

Multiphysics Object Oriented Simulation Environment
https://www.mooseframework.org
GNU Lesser General Public License v2.1
1.6k stars 1k forks source link

Fix TabulatedFluidProperties for multiphase class PorousFlowWaterVapor #27506

Open aikubo opened 2 weeks ago

aikubo commented 2 weeks ago

closes #27457

aikubo commented 2 weeks ago

This is my first moose pull request so I expect there to be some formatting or other errors. thanks for your patience @GiudGiud

Also, I found that using the VScode git interface to try and push and pull causes errors. So I recommend other newbies do it all on terminal.

GiudGiud commented 2 weeks ago

so there's these two formatting errors which should be fairly easy to fix

##########################################################################
ERROR: Your patch does not contain a valid ticket reference! (i.e. #1234)
Merge branch 'tfp_water_vapor' of https://github.com/aikubo/moose into test
style patch
remove debug statements
add tests in PorousFlow fluidstate
Add functions so TabulatedFluidProperties works with multiphase class PorousFlowWaterVapor

ERROR: The following files contain trailing whitespace after applying your patch:
    modules/porous_flow/test/tests/fluidstate/water_vapor_tab.i

Run the "delete_trailing_whitespace.sh" script in your $MOOSE_DIR/scripts directory.
##########################################################################
GiudGiud commented 2 weeks ago

let's remove the merge commit. We prefer rebase. You can remove it by:

git reset --hard a6a0e70
git fetch upstream
git rebase upstream/devel
git push origin HEAD:tfp_water_vapor -f

you have more formatting issues

##########################################################################
ERROR: The following files contain trailing whitespace after applying your patch:
    modules/porous_flow/test/tests/fluidstate/water_vapor_tab.i

Run the "delete_trailing_whitespace.sh" script in your $MOOSE_DIR/scripts directory.
##########################################################################

if you are using VSCode you need to use our user settings, and it wont be a problem ever again https://mooseframework.inl.gov/help/development/VSCode.html

aikubo commented 2 weeks ago

let's remove the merge commit. We prefer rebase.

Thanks. I've fixed this.

if you are using VSCode you need to use our user settings, and it wont be a problem ever again

added these settings

moosebuild commented 2 weeks ago

Job Documentation on 09c8b1d wanted to post the following:

View the site here

This comment will be updated on new commits.

GiudGiud commented 1 week ago

cleaned up the commit history a little. let's see if tests pass

aikubo commented 1 week ago

@cpgr @GiudGiud

I went back to my original input file from #27159. And I get a different error using tabulated fluid properties than from the test input file I made water_vapor_tab.i. And I'm confused...

*** ERROR ***
/home/akh/myprojects/moose_projects/dikes/sims/multi/MVE_s_from_h_p_error.i:186.5:
The following error occurred in the FluidProperties 'water97' of type Water97FluidProperties.

The fluid properties class 'Water97FluidProperties' has not implemented the method below, which computes derivatives of fluid properties with regards to the flow variables. If your application requires this method, you must either implement it or use a different fluid properties  class.

virtual void SinglePhaseFluidProperties::s_from_h_p(libMesh::Real, libMesh::Real, libMesh::Real&, libMesh::Real&, libMesh::Real&) const

You can avoid this error by neglecting the unimplemented derivatives of fluid properties by setting the 'allow_imperfect_jacobians' parameter

Stack frames: 25
0: libMesh::print_trace(std::ostream&)
1: moose::internal::mooseErrorRaw(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >)
2: MooseBase::callMooseError(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, bool) const
3: /home/akh/myprojects/moose_projects/moose/modules/fluid_properties/lib/libfluid_properties-opt.so.0(+0xbbd84) [0x7f0b716bbd84]
4: /home/akh/myprojects/moose_projects/moose/modules/fluid_properties/lib/libfluid_properties-opt.so.0(+0x13fce4) [0x7f0b7173fce4]
5: /home/akh/myprojects/moose_projects/moose/modules/fluid_properties/lib/libfluid_properties-opt.so.0(+0x143574) [0x7f0b71743574]
6: SinglePhaseFluidProperties::T_from_p_h(double, double, double&, double&, double&) const
7: /home/akh/myprojects/moose_projects/moose/modules/fluid_properties/lib/libfluid_properties-opt.so.0(+0x126f35) [0x7f0b71726f35]
8: PorousFlowWaterVapor::thermophysicalProperties(MetaPhysicL::DualNumber<double, MetaPhysicL::SemiDynamicSparseNumberArray<double, unsigned long, MetaPhysicL::NWrapper<64ul> >, true> const&, MetaPhysicL::DualNumber<double, MetaPhysicL::SemiDynamicSparseNumberArray<double, unsigned long, MetaPhysicL::NWrapper<64ul> >, true> const&, unsigned int, FluidStatePhaseEnum&, std::vector<FluidStateProperties, std::allocator<FluidStateProperties> >&) const
9: PorousFlowWaterVapor::thermophysicalProperties(double, double, unsigned int, std::vector<FluidStateProperties, std::allocator<FluidStateProperties> >&) const
10: PorousFlowFluidStateBaseMaterialTempl<false>::initQpStatefulProperties()
11: PorousFlowFluidStateSingleComponentTempl<false>::initQpStatefulProperties()
12: MaterialBase::initStatefulProperties(unsigned int)
13: MaterialPropertyStorage::initStatefulProps(unsigned int, std::vector<std::shared_ptr<MaterialBase>, std::allocator<std::shared_ptr<MaterialBase> > > const&, unsigned int, libMesh::Elem const&, unsigned int)
14: ThreadedElementLoopBase<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> >::operator()(libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> const&, bool)
15: /home/akh/myprojects/moose_projects/moose/framework/libmoose-opt.so.0(+0x20e16d0) [0x7f0b74ee16d0]
16: FEProblemBase::initElementStatefulProps(libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> const&, bool)
17: FEProblemBase::initialSetup()
18: Transient::init()
19: MooseApp::executeExecutioner()
20: MooseApp::run()
21: main
22: /lib/x86_64-linux-gnu/libc.so.6(+0x29d90) [0x7f0b68c29d90]
23: __libc_start_main
24: ../../dikes-opt(+0x33ef) [0x557abb8ce3ef]
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=1
:
system msg for write_line failure : Bad file descriptor

This is an issue with the initial condition for enthalpy. Basically if you're not in the two phase zone, it gives this error. But I don't know why because I don't think any of these objects ever call s_from_h_p and it ONLY happens when using Tabulated Fluid Properties.

Here's a MVE:

[Mesh]
    [gen]
      type = GeneratedMeshGenerator
      dim = 2

    []
  []

  [GlobalParams]
    PorousFlowDictator = dictator
  []

  [AuxVariables]
    [pressure_gas]
      order = CONSTANT
      family = MONOMIAL
    []
    [pressure_water]
      order = CONSTANT
      family = MONOMIAL
    []
    [enthalpy_gas]
      order = CONSTANT
      family = MONOMIAL
    []
    [enthalpy_water]
      order = CONSTANT
      family = MONOMIAL
    []
    [saturation_gas]
      order = CONSTANT
      family = MONOMIAL
    []
    [saturation_water]
      order = CONSTANT
      family = MONOMIAL
    []
    [density_water]
      order = CONSTANT
      family = MONOMIAL
    []
    [density_gas]
      order = CONSTANT
      family = MONOMIAL
    []
    [viscosity_water]
      order = CONSTANT
      family = MONOMIAL
    []
    [viscosity_gas]
      order = CONSTANT
      family = MONOMIAL
    []
    # [temperature]
    #   order = CONSTANT
    #   family = MONOMIAL
    # []
  []

  [AuxKernels]
    [enthalpy_water]
      type = PorousFlowPropertyAux
      variable = enthalpy_water
      property = enthalpy
      phase = 0
      execute_on = 'initial timestep_end'
    []
    [enthalpy_gas]
      type = PorousFlowPropertyAux
      variable = enthalpy_gas
      property = enthalpy
      phase = 1
      execute_on = 'initial timestep_end'
    []
    [pressure_water]
      type = PorousFlowPropertyAux
      variable = pressure_water
      property = pressure
      phase = 0
      execute_on = 'initial timestep_end'
    []
    [pressure_gas]
      type = PorousFlowPropertyAux
      variable = pressure_gas
      property = pressure
      phase = 1
      execute_on = 'initial timestep_end'
    []
    [saturation_water]
      type = PorousFlowPropertyAux
      variable = saturation_water
      property = saturation
      phase = 0
      execute_on = 'initial timestep_end'
    []
    [saturation_gas]
      type = PorousFlowPropertyAux
      variable = saturation_gas
      property = saturation
      phase = 1
      execute_on = 'initial timestep_end'
    []
    [density_water]
      type = PorousFlowPropertyAux
      variable = density_water
      property = density
      phase = 0
      execute_on = 'initial timestep_end'
    []
    [density_gas]
      type = PorousFlowPropertyAux
      variable = density_gas
      property = density
      phase = 1
      execute_on = 'initial timestep_end'
    []
    [viscosity_water]
      type = PorousFlowPropertyAux
      variable = viscosity_water
      property = viscosity
      phase = 0
      execute_on = 'initial timestep_end'
    []
    [viscosity_gas]
      type = PorousFlowPropertyAux
      variable = viscosity_gas
      property = viscosity
      phase = 1
      execute_on = 'initial timestep_end'
    []
  []

  [UserObjects]
    [dictator]
      type = PorousFlowDictator
      porous_flow_vars = 'pliquid h'
      number_fluid_phases = 2
      number_fluid_components = 1
    []
    [pc] #from porousflow/test/tests/fluidstate/watervapor.i
      type = PorousFlowCapillaryPressureBC
      pe = 1e5
      lambda = 2
      pc_max = 1e6
    []
    [fs]
      type = PorousFlowWaterVapor
      water_fp = water_tab
      capillary_pressure = pc
    []
  []

  [Variables]
    [pliquid]
      initial_condition = 1e6
    []
    [h]
      initial_condition = 8e4 # if you change it to 8e5 like in water_vapor_tab.i it works fine
    []
  []

  [Kernels]
    [mass]
      type = PorousFlowMassTimeDerivative
      variable = pliquid
    []
    [heat]
      type = PorousFlowEnergyTimeDerivative
      variable = h

    []
  []

  [FluidProperties]
    [water97]
      type = Water97FluidProperties    # IAPWS-IF97
    []
    [water_tab]
      type = TabulatedBicubicFluidProperties
      fp = water97
      fluid_property_file = 'water_extended.csv'
      p_h_variables = true
    []

  []

  [Materials]
      [watervapor]
        type = PorousFlowFluidStateSingleComponent
        porepressure = pliquid
        enthalpy = h
        capillary_pressure = pc
        fluid_state = fs
        temperature_unit = Kelvin
      []

      [porosity_wallrock]
        type = PorousFlowPorosityConst
        porosity = 0.1
      []
      [permeability]
        type = PorousFlowPermeabilityConst
        permeability = '1E-13 0 0   0 1E-13 0   0 0 1E-13'
      []
      [relperm_water] # from watervapor.i
        type = PorousFlowRelativePermeabilityCorey
        n = 2
        phase = 0
      []
      [relperm_gas]  # from watervapor.i
        type = PorousFlowRelativePermeabilityCorey
        n = 3
        phase = 1
      []
      [internal_energy]
        type = PorousFlowMatrixInternalEnergy
        density = 2500 # kg/m^3
        specific_heat_capacity = 1200 # J/kg/K
      []

  []

  [Preconditioning]
    [smp]
      type = SMP
      full = true
      petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
      petsc_options_value = ' lu       mumps'
    []
  []

  [Executioner]
    type = Transient
    solve_type = NEWTON
    end_time = 100
  []
GiudGiud commented 1 week ago

so you forwarded s_from_h_p to Water97, but is it defined there? Is it defined as s_from_p_h by any chance?

aikubo commented 1 week ago

so you forwarded s_from_h_p to Water97, but is it defined there?

It isn't defined in Water97 and I can add it but I find it odd because I'm not sure why this error occurs.

Is it defined as s_from_p_h by any chance?

no it is not but it was a good idea to check that

GiudGiud commented 1 week ago

the backtrace explains the story: -something in porousflow water calls T_from_p_h

so two solutions:

aikubo commented 1 week ago

Could you just implement it as basically a copy of s_from_p_T? Or is more complexity necessary?

template <typename T>
T
Water97FluidProperties::s_from_h_p_template(const T & enthalpy, const T & pressure) const
{
  T temperature = T_from_h_p(enthalpy, pressure);
  return s_from_p_T_template(pressure, temperature);
}

template <typename T>
void
Water97FluidProperties::s_from_h_p_template(
    const T & enthalpy, const T & pressure, T & s, T & ds_dh, T & ds_dp) const
{
  auto functor = [this](const auto & enthalpy, const auto & pressure)
  { return this->s_from_p_T_template(enthalpy, pressure); };

  xyDerivatives(enthalpy, pressure, s, ds_dh, ds_dp, functor);
}

And another question, looking closely there's T_from_p_h and T_from_h_p in Water97FluidProperties (one of which I added)
And in SinglePhaseFluidProperties there's both. But only T_from_h_p in Tabulated Fluid Properties. Aren't these the same? Which should it be?

GiudGiud commented 1 week ago 
template <typename T>
T
Water97FluidProperties::s_from_h_p_template(const T & enthalpy, const T & pressure) const
{
  T temperature = T_from_h_p(enthalpy, pressure);
  return s_from_p_T_template(pressure, temperature);
}

this would be fine. As long as T_from_h_p is implemented. Seems like it is.

And another question, looking closely there's T_from_p_h and T_from_h_p in Water97FluidProperties (one of which I added)

Yes this is unfortunate. We only want T_from_p_h, but they both got implemented. See in SPFP.h:

  propfunc(T, h, p)  // temporary, until uniformization
  propfuncWithDefault(T, p, h)

T_from_p_h is the official one

GiudGiud commented 1 week ago

but if T_from_h_p is implemented, T_from_p_h is definitely the faster and easier to implement

aikubo commented 1 week ago

I've added s_from_h_p and it works fine but I kept getting issues with 

*** ERROR ***
/home/akh/myprojects/moose_projects/dikes/MVE_s_from_h_p_error.i:183.5:
The following error occurred in the FluidProperties 'water_tab' of type TabulatedBicubicFluidProperties.

You must construct pT from vh tables when calling e_from_v_h.

When you turn both construct_pT_from_ve=true and construct_pT_from_vh=true everything works fine. This only occurs in certain regions of the phase diagram for p,h variables in porousflow. I think it might make sense to turn on those flags automatically as part of p_h_variables flag in TabulatedFluidProperties. Or if not to note it somewhere.

GiudGiud commented 1 week ago

. I think it might make sense to turn on those flags automatically as part of p_h_variables flag in TabulatedFluidProperties.

When does this flag get set? By the user? If so then sure.

Or if not to note it somewhere.

That's fine too. I'm guessing that you are hitting these 'new' calls because:

aikubo commented 1 day ago

I was having a hard time tracking down some of the issues I was having and I'll do my best to summarize

I tried a input file in the Liquid state and it wasn't running at all (MVE below). I added the s_from_h_p to Water97FluidProperties and accordingly to TabulatedFluidProperties. Then the issue was that TabulatedFluidProperties was calling some of the v,h and v,e functions so it needed to construct that table. So I added it in that p_h_variables automatically turns those on. p_h_variables is set by the user in the input file.

After that the input file ran but the temperature was wrong. It gave a temperature of 299 K when it should be 377 K. I tracked this to PorousFlowWaterVapor calling T_from_p_h but the AD version which wasn't implemented in TabulatedFluidProperties, only the Real version.

I also had to add a ADReal T_from_p_h version to WaterFluidProperties because the current version is FPADReal T_from_p_h_ad.

The way I implemented ADReal T_from_p_h is imho terrible but I had a hard time compiling between the FPADReal and ADReal types. It's unclear to me when to use FPADReal, it seems to be used mostly for internal FP code. And as far as I can tell the T_from_p_h function is only directly called outside the FluidProperties module in PorousFlowWaterVapor . What is the correct usage and correct way to convert between types?

Finally, I cleaned up the commit history.

Thanks for your patience bearing with a MOOSE and C++ novice.

Here's the liquid state MVE which works now:

# at 1MPa (1e6 Pa) and 440 kJ/kg (440000) T should be 377 K

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[AuxVariables]
  [pressure_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [pressure_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [temperature]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [enthalpy_water]
    type = PorousFlowPropertyAux
    variable = enthalpy_water
    property = enthalpy
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [enthalpy_gas]
    type = PorousFlowPropertyAux
    variable = enthalpy_gas
    property = enthalpy
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_water]
    type = PorousFlowPropertyAux
    variable = pressure_water
    property = pressure
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [pressure_gas]
    type = PorousFlowPropertyAux
    variable = pressure_gas
    property = pressure
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [density_water]
    type = PorousFlowPropertyAux
    variable = density_water
    property = density
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = PorousFlowPropertyAux
    variable = density_gas
    property = density
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [viscosity_water]
    type = PorousFlowPropertyAux
    variable = viscosity_water
    property = viscosity
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [viscosity_gas]
    type = PorousFlowPropertyAux
    variable = viscosity_gas
    property = viscosity
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [temperature]
    type = PorousFlowPropertyAux
    variable = temperature
    property = temperature
    execute_on = 'initial timestep_end'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pliquid h'
    number_fluid_phases = 2
    number_fluid_components = 1
  []
  [pc] #from porousflow/test/tests/fluidstate/watervapor.i
    type = PorousFlowCapillaryPressureBC
    pe = 1e5
    lambda = 2
    pc_max = 1e6
  []
  [fs]
    type = PorousFlowWaterVapor
    water_fp = water_tab
    capillary_pressure = pc
  []
[]

[Variables]
  [pliquid]
    initial_condition = 1e6
  []
  [h]
    initial_condition = 440e3
  []
[]

[Kernels]
  [mass]
    type = PorousFlowMassTimeDerivative
    variable = pliquid
  []
  [heat]
    type = PorousFlowEnergyTimeDerivative
    variable = h
  []
[]

[FluidProperties]
  [water97]
    type = Water97FluidProperties # IAPWS-IF97
  []
  [water_tab]
    type = TabulatedBicubicFluidProperties
    fp = water97
    p_h_variables = true
  []
[]

[Materials]
  [watervapor]
    type = PorousFlowFluidStateSingleComponent
    porepressure = pliquid
    enthalpy = h
    capillary_pressure = pc
    fluid_state = fs
    temperature_unit = Kelvin
  []

  [porosity_wallrock]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-13 0 0   0 1E-13 0   0 0 1E-13'
  []
  [relperm_water] # from watervapor.i
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm_gas] # from watervapor.i
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500 # kg/m^3
    specific_heat_capacity = 1200 # J/kg/K
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Postprocessors]
  [density_water]
    type = ElementAverageValue
    variable = density_water
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = ElementAverageValue
    variable = density_gas
    execute_on = 'initial timestep_end'
  []
  [viscosity_water]
    type = ElementAverageValue
    variable = viscosity_water
    execute_on = 'initial timestep_end'
  []
  [viscosity_gas]
    type = ElementAverageValue
    variable = viscosity_gas
    execute_on = 'initial timestep_end'
  []
  [enthalpy_water]
    type = ElementAverageValue
    variable = enthalpy_water
    execute_on = 'initial timestep_end'
  []
  [enthalpy_gas]
    type = ElementAverageValue
    variable = enthalpy_gas
    execute_on = 'initial timestep_end'
  []
  [sg]
    type = ElementAverageValue
    variable = saturation_gas
    execute_on = 'initial timestep_end'
  []
  [sw]
    type = ElementAverageValue
    variable = saturation_water
    execute_on = 'initial timestep_end'
  []
  [pwater]
    type = ElementAverageValue
    variable = pressure_water
    execute_on = 'initial timestep_end'
  []
  [pgas]
    type = ElementAverageValue
    variable = pressure_gas
    execute_on = 'initial timestep_end'
  []
  [temperature]
    type = ElementAverageValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [enthalpy]
    type = ElementAverageValue
    variable = h
    execute_on = 'initial timestep_end'
  []
  [liquid_mass]
    type = PorousFlowFluidMass
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [vapor_mass]
    type = PorousFlowFluidMass
    phase = 1
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  file_base = water_vapor_twophase_tab
  csv = true
[]
GiudGiud commented 16 hours ago

Sounds good. I ll check this out next week. Then I will be on travel so I may not be able to help this through. If you mind the delay we will have to bring in another reviewer.

I don't think we need to use the FP version of the ADReal. I get by with using ADReal everywhere