Divergence of phi not available after reinit

grmnptr commented 3 months ago

Bug Description

Even though some kernels request the divergence of vector basis functions (the one we saw issues with wereNEDELEC_ONE and LAGRANGE_VEC), it is not available after reinit and in opt mode we hit a segfault while in debug mode we hit an assertion with calculations_started as true and calculate_div_phi as false at the top of the stack below:

thread #1, queue = 'com.apple.main-thread', stop reason = breakpoint 1.1 2.1 3.1 4.1
  * frame #0: 0x0000000107ec9444 libmpi.12.dylib`MPI_Abort
    frame #1: 0x00000001373157d0 libmesh_dbg.0.dylib`libMesh::libmesh_terminate_handler() + 320
    frame #2: 0x000000018a3990cc libc++abi.dylib`std::__terminate(void (*)()) + 16
    frame #3: 0x000000018a39c348 libc++abi.dylib`__cxxabiv1::failed_throw(__cxxabiv1::__cxa_exception*) + 88
    frame #4: 0x000000018a39c28c libc++abi.dylib`__cxa_throw + 308
    frame #5: 0x0000000110c9c6b4 libmoose-dbg.0.dylib`libMesh::FEGenericBase<libMesh::VectorValue<double>>::get_div_phi(this=0x00000001001486d0) const at fe_base.h:262:5
    frame #6: 0x0000000110ca2e64 libmoose-dbg.0.dylib`Assembly::reinitFE(this=0x000000013389e400, elem=0x0000600001739b00) at Assembly.C:800:81
    frame #7: 0x0000000110d14798 libmoose-dbg.0.dylib`Assembly::reinit(this=0x000000013389e400, elem=0x0000600001739b00) at Assembly.C:1820:3
    frame #8: 0x000000010f908a18 libmoose-dbg.0.dylib`FEProblemBase::prepare(this=0x000000013389a800, elem=0x0000600001739b00, tid=0) at FEProblemBase.C:1491:24
    frame #9: 0x0000000110671ea0 libmoose-dbg.0.dylib`ComputeInitialConditionThread::operator()(this=0x000000016fdf5e88, range=0x00006000002686e0) at ComputeInitialConditionThread.C:75:17
    frame #10: 0x000000010f8bfa4c libmoose-dbg.0.dylib`void libMesh::Threads::parallel_reduce<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*>, ComputeInitialConditionThread>(range=0x00006000002686e0, body=0x000000016fdf5e88) at threads_pthread.h:370:5
    frame #11: 0x000000010f8bf2a4 libmoose-dbg.0.dylib`FEProblemBase::projectSolution(this=0x000000013389a800) at FEProblemBase.C:3254:3
    frame #12: 0x000000010f9227d4 libmoose-dbg.0.dylib`FEProblemBase::initialSetup(this=0x000000013389a800) at FEProblemBase.C:961:5
    frame #13: 0x000000010ebcf0e0 libmoose-dbg.0.dylib`Transient::init(this=0x0000000133894a00) at Transient.C:225:12
    frame #14: 0x0000000110ce714c libmoose-dbg.0.dylib`MooseApp::executeExecutioner(this=0x0000000133889618) at MooseApp.C:1157:19
    frame #15: 0x0000000110ce1300 libmoose-dbg.0.dylib`MooseApp::run(this=0x0000000133889618) at MooseApp.C:1541:5
    frame #16: 0x000000010001753c test_kernels-dbg`void Moose::main<test_kernelsTestApp>(argc=3, argv=0x000000016fdfdc80) at MooseMain.h:47:8
    frame #17: 0x00000001000174bc test_kernels-dbg`main(argc=3, argv=0x000000016fdfdc80) at main.C:17:3
    frame #18: 0x000000018a05a0e0 dyld`start + 2360

which interrupts here:

if (_need_div.find(fe_type) != _need_div.end())
      fesd._div_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe.get_div_phi()));

Based on a discussion with @roystgnr this could be fixed by adding the conditions in the prerequests too.

Steps to Reproduce

Run the following input files:

angular_frequency = 500
electric_conductivity = 0.01

[Mesh]
    [fmg]
        any mesh
    []
[]

[Variables]
    [A_real]
        family = LAGRANGE_VEC
    []
    [V_real]
        family = LAGRANGE
    []

[]

[Functions]
    [BC_V_left_port]
        type = ParsedFunction
        expression = '10.0 * cos(2*pi*500*t)'
    []
    [BC_V_right_port]
        type = ParsedFunction
        expression = '0'
    []
[]

[Kernels]  
    ## Coupled kernels for V_real & V_imag
    [Coupled_real_V_real]
        type = DivField
        variable = V_real
        coupled_vector_variable = A_real
    []
    [Time_deriv_V_real]
        type = TimeDerivative
        variable = V_real
    []
    [Diffusion_A_real]
        type = VectorDiffusion
    variable = A_real
    []
[]

[BCs]
    [Left_port_V_real]
        type = FunctionDirichletBC
        variable = V_real
        function = BC_V_left_port
        boundary = left
    []
    [Right_port_V_real]
        type = FunctionDirichletBC
        variable = V_real
        function = BC_V_right_port
        boundary = right
    []
[]

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

[Executioner]
    type = Transient
    solve_type = 'NEWTON'
    nl_abs_tol = 1e-8
    nl_rel_tol = 1e-8
    nl_max_its = 100
    dt = 0.001
    num_steps = 1
    petsc_options_iname = '-pc_type'
    petsc_options_value = 'lu'
    automatic_scaling = true
    compute_scaling_once = false
    scaling_group_variables = 'V_real V_imag'
    [TimeIntegrator]
        type = NewmarkBeta
    []
[]

[Outputs]
    exodus = true
    vtk = true 
    csv = True
    print_linear_residuals = true
[]

[Debug]
    show_var_residual_norms = false
[]

Impact

Code crashes.

[Optional] Diagnostics

No response

roystgnr commented 3 months ago

Hey, are you guys using an up-to-date MOOSE here? I've got a fix for a related bug with get_curl_phi(), along with a nice refactoring, so I'm glad I looked at this, but I was actually only able to replicate the get_div_phi() issue because the first MOOSE build I grabbed was old enough to predate #26477

lindsayad commented 3 months ago

Small kudos to me I guess 😆

grmnptr commented 3 months ago

Nono, large kudos to you! Thanks @roystgnr too for looking into this!

lindsayad commented 3 months ago

Were you guys using an up-to-date MOOSE per @roystgnr's question?

grmnptr commented 3 months ago

I don't think so. I told the intern (not sure who is the mentor) to update his MOOSE and let us know if he sees more issues. I think we can probably close this.

lindsayad commented 3 months ago

We can keep it open a bit longer to see if we hear back

idaholab / moose