Differences between convection integrators

[sbugenha]

Hi All,

I'm implementing the incompressible Navier Stokes equations using the IDA solver from Sundials (see PR #3690) and am confused about which convection integrator I should be using in MFEM. There are 3 options - 2 nonlinear form integrators (VectorConvectionNLFIntegrator and ConvectiveVectorConvectionNLFIntegrator) and the bilinear form ConvectionIntegrator. Applying the bilinear form version using:

    k = new BilinearForm(spaces[2]);
    k->AddDomainIntegrator(new ConvectionIntegrator(w_coeff, rho));
    k->Assemble();
    k->Finalize();
    K = k->LoseMat();
    Kmat = new BlockMatrix(componentOffsets);
    Kmat->SetBlock(0, 0, K);
    Kmat->SetBlock(1, 1, K);
    Kmonolithic = Kmat->CreateMonolithic();

gives me similar results to applying the nonlinear form using:

    N = new NonlinearForm(spaces[0]);
    N->AddDomainIntegrator(new VectorConvectionNLFIntegrator(rho_coeff));
    //N->AddDomainIntegrator(new ConvectiveVectorConvectionNLFIntegrator(rho_coeff));
        gradN = dynamic_cast<SparseMatrix*>(&N->GetGradient(u));

where Kmonolithic and gradN are the gradient matrices of the bilinear and nonlinear forms, respectively. Here, spaces[0] is a 2D vector space whereas spaces[2] is a 2D scalar space. However, applying the ConvectiveVectorConvectionNLFIntegrator in the nonlinear form (commented line above) gives me very different results from the ConvectionIntegrator and VectorConvectionNLFIntegrator.

What are the differences between these approaches? What approach is most appropriate for IDA? IDA uses a modified Newton iteration where the Jacobian is fixed and usually out of date (and only updated as frequently as necessary) throughout the nonlinear iterations.

Thank you for your help.

-Scott

[jandrej]

I can only vouch for the VectorConvectionNLFIntegrator and its capabilities to compute the gradient for a Newton method as I have not used ConvectiveVectorConvectionNLFIntegrator.

ConvectionIntegrator is not the right integrator to use for Navier-Stokes as it won't capture off-diagonal terms. Depending on your velocity field this may line up with VectorConvectionNLFIntegrator.

Using IDA for Navier-Stokes can work but the hard part is not the time integration itself but the solution of the linear system. Specifically the linearization of VectorConvectionNLFIntegrator which appears in A(1,1) of your linearization during a Newton iteration in combination with Reynolds numbers > 100.

Keep us posted.

[sbugenha]

Thanks @jandrej. Just to clarify, the action of the ConvectionIntegrator and VectorConvectionNLFIntegrator obtained by Mult should be the same, while the gradient matrices produced by the integrators differs, is this correct? I.e. there should be no difference in the residual evaluation, just the Jacobian? Also, if the convection term is linearized by lagging the convection coefficient (e.g. with the Oseen iteration), then would the ConvectionIntegrator be appropriate?

Using IDA for Navier-Stokes can work but the hard part is not the time integration itself but the solution of the linear system.

Totally agree. I manage to get reasonable iteration number solving the linear system using the Pressure-Convection-Diffusion block preconditioner which is relatively independent of mesh and Reynolds number when I use direct solves on the momentum and pressure-Poisson blocks, but I've been having some issues with convergence solving the momentum block using GMRES with BoomerAMG preconditioner. Part of it may be that I was using the ConvectiveVectorConvectionNLFIntegrator rather than the VectorConvectionNLFIntegrator. But I also suspect stability issues with convection dominated flow may be to blame and am working on implementing SUPG stabilization to see if this helps.

I'll definitely keep you posted. Implementing the stabilization doesn't seem to be as straightforward as I was hoping, so I expect I'll need to reach out to you for help again as I work through this.

Thanks, Scott

[jandrej]

the action of the ConvectionIntegrator and VectorConvectionNLFIntegrator obtained by Mult should be the same, while the gradient matrices produced by the integrators differs, is this correct?

Yes. Making it work involves more work though. You need component wise operations.

Also, if the convection term is linearized by lagging the convection coefficient (e.g. with the Oseen iteration), then would the ConvectionIntegrator be appropriate?

Yes, but see above.

but I've been having some issues with convergence solving the momentum block using GMRES with BoomerAMG preconditioner

This is expected. The larger the Reynolds number, the more the operator drives into convective regime.

Implementing the stabilization doesn't seem to be as straightforward as I was hoping

If you're looking into SUPG/PSPG, you have to dive into the integrators and implement the terms right there.

[stale[bot]]

:warning: This issue or PR has been automatically marked as stale because it has not had any activity in the last month. If no activity occurs in the next week, it will be automatically closed. Thank you for your contributions.