Got an error when using adaptivity

[doquang]

Check these boxes if you have followed the posting rules.

• Q&A Navier Stokes is the most appropriate category 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 (avoid screenshots if possible, triple back quotes before/after pasted text, etc)

Question

When I tried to extend the gallium_melting example to use adaptive remeshing, I encountered an error after a few steps.
Could someone help!

You must call one or more of the FE accessors (e.g. get_phi(), get_dphi(), get_nothing()) _before_ calling reinit()!
Stack frames: 20
0: 0   libmesh_opt.0.dylib                 0x000000010accc9c4 libMesh::print_trace(std::__1::basic_ostream<char, std::__1::char_traits<char>>&) + 840
1: 1   libmesh_opt.0.dylib                 0x000000010acc981c libMesh::MacroFunctions::report_error(char const*, int, char const*, char const*, std::__1::basic_ostream<char, std::__1::char_traits<char>>&) + 296
2: 2   libmesh_opt.0.dylib                 0x000000010adeedac libMesh::FEGenericBase<double>::determine_calculations() + 344
3: 3   libmesh_opt.0.dylib                 0x000000010ad45604 libMesh::FE<2u, (libMesh::FEFamily)31>::reinit(libMesh::Elem const*, std::__1::vector<libMesh::Point, std::__1::allocator<libMesh::Point>> const*, std::__1::vector<double, std::__1::allocator<double>> const*) + 64
4: 4   libmesh_opt.0.dylib                 0x000000010b4445a4 libMesh::FEMContext::elem_fe_reinit(std::__1::vector<libMesh::Point, std::__1::allocator<libMesh::Point>> const*) + 156
5: 5   libmesh_opt.0.dylib                 0x000000010b49684c libMesh::GenericProjector<libMesh::OldSolutionValue<double, &void libMesh::FEMContext::point_value<double>(unsigned int, libMesh::Point const&, double&, double) const>, libMesh::OldSolutionValue<libMesh::VectorValue<double>, &void libMesh::FEMContext::point_gradient<libMesh::VectorValue<double>>(unsigned int, libMesh::Point const&, double&, double) const>, double, libMesh::VectorSetAction<double>>::ProjectInteriors::operator()(libMesh::StoredRange<std::__1::__wrap_iter<libMesh::Elem const* const*>, libMesh::Elem const*> const&) + 1176
6: 6   libmesh_opt.0.dylib                 0x000000010b488410 void libMesh::Threads::parallel_reduce<libMesh::StoredRange<std::__1::__wrap_iter<libMesh::Elem const* const*>, libMesh::Elem const*>, libMesh::GenericProjector<libMesh::OldSolutionValue<double, &void libMesh::FEMContext::point_value<double>(unsigned int, libMesh::Point const&, double&, double) const>, libMesh::OldSolutionValue<libMesh::VectorValue<double>, &void libMesh::FEMContext::point_gradient<libMesh::VectorValue<double>>(unsigned int, libMesh::Point const&, double&, double) const>, double, libMesh::VectorSetAction<double>>::ProjectInteriors>(double const&, libMesh::GenericProjector<libMesh::OldSolutionValue<double, &void libMesh::FEMContext::point_value<double>(unsigned int, libMesh::Point const&, double&, double) const>, libMesh::OldSolutionValue<libMesh::VectorValue<double>, &void libMesh::FEMContext::point_gradient<libMesh::VectorValue<double>>(unsigned int, libMesh::Point const&, double&, double) const>, double, libMesh::VectorSetAction<double>>::ProjectInteriors&, unsigned int) + 112
7: 7   libmesh_opt.0.dylib                 0x000000010b48401c libMesh::GenericProjector<libMesh::OldSolutionValue<double, &void libMesh::FEMContext::point_value<double>(unsigned int, libMesh::Point const&, double&, double) const>, libMesh::OldSolutionValue<libMesh::VectorValue<double>, &void libMesh::FEMContext::point_gradient<libMesh::VectorValue<double>>(unsigned int, libMesh::Point const&, double&, double) const>, double, libMesh::VectorSetAction<double>>::project(libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> const&) + 2016
8: 8   libmesh_opt.0.dylib                 0x000000010b48240c libMesh::System::project_vector(libMesh::NumericVector<double> const&, libMesh::NumericVector<double>&, int, std::__1::optional<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*>>, std::__1::optional<std::__1::vector<unsigned int, std::__1::allocator<unsigned int>>>) const + 2308
9: 9   libmesh_opt.0.dylib                 0x000000010b4819e8 libMesh::System::project_vector(libMesh::NumericVector<double>&, int, std::__1::optional<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*>>, std::__1::optional<std::__1::vector<unsigned int, std::__1::allocator<unsigned int>>>) const + 168
10: 10  libmesh_opt.0.dylib                 0x000000010b45e80c libMesh::System::restrict_vectors() + 368
11: 11  libmesh_opt.0.dylib                 0x000000010b42e91c libMesh::EquationSystems::reinit_solutions() + 1784
12: 12  libmoose-opt.0.dylib                0x00000001061dd7ec FEProblemBase::meshChanged(bool, bool, bool) + 320
13: 13  libmoose-opt.0.dylib                0x00000001061eca9c FEProblemBase::adaptMesh() + 340
14: 14  libmoose-opt.0.dylib                0x0000000105b72198 TransientBase::incrementStepOrReject() + 96
15: 15  libmoose-opt.0.dylib                0x0000000105b70678 TransientBase::execute() + 912
16: 16  libmoose-opt.0.dylib                0x0000000106cb96f8 MooseApp::executeExecutioner() + 368
17: 17  libmoose-opt.0.dylib                0x0000000106cb4f7c MooseApp::run() + 1760
18: 18  tcmoose-opt                         0x0000000102c61908 main + 608
19: 19  dyld                                0x0000000189033e00 start + 6992
[0] ../src/fe/fe_base.C, line 932, compiled May 13 2026 at 12:06:03

##########################################################
# Simulation of Galium Melting Experiment
# Ref: Gau, C., & Viskanta, R. (1986). Melting and solidification of a pure metal on a vertical wall.
# Key physics: melting/solidification, convective heat transfer, natural convection
##########################################################

mu = 1.81e-3
rho_solid = 6093
rho_liquid = 6093
k_solid = 32
k_liquid = 32
cp_solid = 381.5
cp_liquid = 381.5
L = 80160
alpha_b = 1.2e-4
T_solidus = 302.93
T_liquidus = '${fparse T_solidus + 0.1}'
advected_interp_method = 'upwind'
velocity_interp_method = 'rc'
T_cold = 301.15
T_hot = 311.15
Nx = 20
Ny = 10

[GlobalParams]
  rhie_chow_user_object = 'rc'
[]

[UserObjects]
  [rc]
    type = INSFVRhieChowInterpolator
    u = vel_x
    v = vel_y
    pressure = pressure
  []
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 88.9e-3
    ymin = 0
    ymax = 63.5e-3
    nx = ${Nx}
    ny = ${Ny}
  []
[]

[AuxVariables]
  [U]
    type = MooseVariableFVReal
  []
  [fl]
    type = MooseVariableFVReal
    initial_condition = 0.0
  []
  [density]
    type = MooseVariableFVReal
  []
  [th_cond]
    type = MooseVariableFVReal
  []
  [cp_var]
    type = MooseVariableFVReal
  []
  [darcy_coef]
    type = MooseVariableFVReal
  []
  [fch_coef]
    type = MooseVariableFVReal
  []
[]

[AuxKernels]
  [mag]
    type = VectorMagnitudeAux
    variable = U
    x = vel_x
    y = vel_y
  []
  [compute_fl]
    type = NSLiquidFractionAux
    variable = fl
    temperature = T
    T_liquidus = '${T_liquidus}'
    T_solidus = '${T_solidus}'
    execute_on = 'TIMESTEP_END'
  []
  [rho_out]
    type = FunctorAux
    functor = 'rho_mixture'
    variable = 'density'
  []
  [th_cond_out]
    type = FunctorAux
    functor = 'k_mixture'
    variable = 'th_cond'
  []
  [cp_out]
    type = FunctorAux
    functor = 'cp_mixture'
    variable = 'cp_var'
  []
  [darcy_out]
    type = FunctorAux
    functor = 'Darcy_coefficient'
    variable = 'darcy_coef'
  []
  [fch_out]
    type = FunctorAux
    functor = 'Forchheimer_coefficient'
    variable = 'fch_coef'
  []
[]

[Variables]
  [vel_x]
    type = INSFVVelocityVariable
    initial_condition = 0.0
  []
  [vel_y]
    type = INSFVVelocityVariable
    initial_condition = 0.0
  []
  [pressure]
    type = INSFVPressureVariable
  []
  [lambda]
    family = SCALAR
    order = FIRST
  []
  [T]
    type = INSFVEnergyVariable
    initial_condition = '${T_cold}'
    scaling = 1e-4
  []
[]

[FVKernels]
  [mass]
    type = INSFVMassAdvection
    variable = pressure
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = rho_mixture
  []
  [mean_zero_pressure]
    type = FVIntegralValueConstraint
    variable = pressure
    lambda = lambda
    phi0 = 0.0
  []

  [u_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_x
    rho = rho_mixture
    momentum_component = 'x'
  []
  [u_advection]
    type = INSFVMomentumAdvection
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = rho_mixture
    momentum_component = 'x'
  []
  [u_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_x
    mu = ${mu}
    momentum_component = 'x'
  []
  [u_pressure]
    type = INSFVMomentumPressure
    variable = vel_x
    momentum_component = 'x'
    pressure = pressure
  []
  [u_friction]
    type = PINSFVMomentumFriction
    variable = vel_x
    momentum_component = 'x'
    u = vel_x
    v = vel_y
    Darcy_name = 'Darcy_coeff'
    Forchheimer_name = 'Forchheimer_coeff'
    rho = ${rho_liquid}
    mu = ${mu}
    standard_friction_formulation = false
  []
  [u_buoyancy]
    type = INSFVMomentumBoussinesq
    variable = vel_x
    T_fluid = T
    gravity = '0 -9.81 0'
    rho = '${rho_liquid}'
    ref_temperature = ${T_cold}
    momentum_component = 'x'
  []
  [u_gravity]
    type = INSFVMomentumGravity
    variable = vel_x
    gravity = '0 -9.81 0'
    rho = '${rho_liquid}'
    momentum_component = 'x'
  []

  [v_time]
    type = INSFVMomentumTimeDerivative
    variable = vel_y
    rho = rho_mixture
    momentum_component = 'y'
  []
  [v_advection]
    type = INSFVMomentumAdvection
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = rho_mixture
    momentum_component = 'y'
  []
  [v_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_y
    mu = ${mu}
    momentum_component = 'y'
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = vel_y
    momentum_component = 'y'
    pressure = pressure
  []
  [v_friction]
    type = PINSFVMomentumFriction
    variable = vel_y
    momentum_component = 'y'
    u = vel_x
    v = vel_y
    Darcy_name = 'Darcy_coeff'
    Forchheimer_name = 'Forchheimer_coeff'
    rho = ${rho_liquid}
    mu = ${mu}
    standard_friction_formulation = false
  []
  [v_buoyancy]
    type = INSFVMomentumBoussinesq
    variable = vel_y
    T_fluid = T
    gravity = '0 -9.81 0'
    rho = '${rho_liquid}'
    ref_temperature = ${T_cold}
    momentum_component = 'y'
  []
  [v_gravity]
    type = INSFVMomentumGravity
    variable = vel_y
    gravity = '0 -9.81 0'
    rho = '${rho_liquid}'
    momentum_component = 'y'
  []

  [T_time]
    type = INSFVEnergyTimeDerivative
    variable = T
    rho = rho_mixture
    dh_dt = dh_dt
  []
  [energy_advection]
    type = INSFVEnergyAdvection
    variable = T
    velocity_interp_method = ${velocity_interp_method}
    advected_interp_method = ${advected_interp_method}
  []
  [energy_diffusion]
    type = FVDiffusion
    coeff = k_mixture
    variable = T
  []
  [energy_source]
    type = NSFVPhaseChangeSource
    variable = T
    L = ${L}
    liquid_fraction = fl
    T_liquidus = ${T_liquidus}
    T_solidus = ${T_solidus}
    rho = 'rho_mixture'
  []

[]

[FVBCs]
  [walls-u]
    type = INSFVNoSlipWallBC
    boundary = 'left right top bottom'
    variable = vel_x
    function = 0
  []
  [walls-v]
    type = INSFVNoSlipWallBC
    boundary = 'left right top bottom'
    variable = vel_y
    function = 0
  []
  [hot_wall]
    type = FVDirichletBC
    variable = T
    value = '${T_hot}'
    boundary = 'left'
  []
  [cold_wall]
    type = FVDirichletBC
    variable = T
    value = '${T_cold}'
    boundary = 'right'
  []
[]

[FunctorMaterials]
  [ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    rho = rho_mixture
    cp = cp_mixture
    temperature = 'T'
  []
  [eff_cp]
    type = NSFVMixtureFunctorMaterial
    phase_2_names = '${cp_solid} ${k_solid} ${rho_solid}'
    phase_1_names = '${cp_liquid} ${k_liquid} ${rho_liquid}'
    prop_names = 'cp_mixture k_mixture rho_mixture'
    phase_1_fraction = fl
  []
  [mushy_zone_resistance]
    type = INSFVMushyPorousFrictionFunctorMaterial
    liquid_fraction = 'fl'
    mu = '${mu}'
    rho_l = '${rho_liquid}'
    dendrite_spacing_scaling = 1e-1
  []
  [friction]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'Darcy_coeff Forchheimer_coeff'
    prop_values = 'darcy_coef darcy_coef darcy_coef fch_coef fch_coef fch_coef'
  []
  [const_functor]
    type = ADGenericFunctorMaterial
    prop_names = 'alpha_b'
    prop_values = '${alpha_b}'
  []
[]

[Adaptivity]
   max_h_level = 1
   marker = temp_maker

  [Markers]
    [temp_maker]
      type = ValueRangeMarker
      lower_bound = ${T_solidus}
      upper_bound = 10000
      variable = T
    []
  []
[]

[Executioner]
  type = Transient

  # Time-stepping parameters
  start_time = 0.0
  end_time = 200.0

  [TimeStepper]
    type = IterationAdaptiveDT
    # Raise time step often but not by as much
    # There's a rough spot for convergence near 10% fluid fraction
    optimal_iterations = 15
    growth_factor = 1.5
    dt = 0.1
  []

  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu NONZERO'
  nl_rel_tol = 1e-6
  nl_max_its = 30
  line_search = 'none'
[]

[Postprocessors]
  [ave_p]
    type = ElementAverageValue
    variable = 'pressure'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [ave_fl]
    type = ElementAverageValue
    variable = 'fl'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [ave_T]
    type = ElementAverageValue
    variable = 'T'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[VectorPostprocessors]
  [vel_x]
    type = ElementValueSampler
    variable = 'vel_x fl'
    sort_by = 'x'
  []
[]

[Outputs]
  exodus = true
  csv = false
[]

Additional information

Mesh size and type:
Reynolds number:
Discretization (finite element CG/DG, finite volume, etc):
Models (turbulence, porous media, etc):
Solver method (fully coupled, segregated, multiapps, etc):
Base input you started from:

[GiudGiud]

Hello

Thanks for reporting the bug

@roystgnr I can reproduce this.
Libmesh is trying to project the scalar variable (lambda) and it's failing because we have not requested everything we'd need for that because we are using finite volume.

Do you prefer:

• building a variable exclusion list in libmesh in system::restrict vector before calling projections?
• use requesting the phi computation for scalar variables in MOOSE

I think the latter might be easier

[roystgnr]

Time step 78 it is, yeah, I can reproduce this too.

we could avoid this from the MOOSE side if we requested more things be computed in FV

I don't think we could. Requests aren't global variables; I'm not a total hack! Adding more requests to an FE object within MOOSE isn't going to add them to the FE objects within the FEMContext objects created by GenericProjector.

Ping me on Monday? I'm stuck on my MacBook right now and it'll be so much nicer to debug back at my workstation.

[GiudGiud]

Le requested reminder @roystgnr

[doquang]

@roystgnr , @GiudGiud - I would like to follow up on this issue. Has there been any progress since our last discussion?

[GiudGiud]

hello @doquang
I made a branch with AI last week that bypasses the issue but I don't think it's quite the right fix
https://github.com/GiudGiud/libmesh/tree/PR_scalar_amr

feel free to use it locally if needed

[GiudGiud]

Opened #33366 so we track this