Hybrid Continuous/Discontinuous Galerkin Finite Element Navier Stokes

In (Pandare and Luo, 2016) the authors describe a scheme in which the velocity is discretized using a reconstructed discontinuous Galerkin (DG) method while the pressure is discretized using continuous Galerkin (CG) finite elements. While MOOSE has the ability to perform rDG(P0P1), otherwise known as a second order accurate finite volume method (see Incompressible Finite Volume Navier Stokes), it does not support arbitrarily high-order reconstructions. It does, on the other hand, support arbitrarily high-order discontinuous Galerkin bases. In that vein we have implemented and tested a hybrid CG-DG scheme. A nice advantage of the hybrid CG-DG scheme is that it is intrinsically LBB stable, allowing equal polynomial order bases for the velocity and pressure (Pandare and Luo, 2016). Additionally, because the velocity is discretized with a DG scheme, the momentum advection flux can be naturally upwinded without a complex upwinding scheme like as Streamline-Upwind Petrov-Galerkin (SUPG) (see Continuous Galerkin Finite Element Navier Stokes).

In the following we run through an example of a hybrid CG-DG input, in which we are using the scheme to solve a lid driven cavity problem with a Reynolds number of 200.

The first block is for the mesh. We create a simple square and add a nodeset which we will use to pin the pressure, removing its nullspace (the discretization only involves the gradient of pressure so without something like a pressure pin, the pressure is only defined up to a constant).

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${l}
    ymin = 0
    ymax = ${l}
    nx = 20
    ny = 20
  []
  [corner_node]
    type = ExtraNodesetGenerator
    new_boundary = 'pinned_node'
    nodes = '0'
    input = gen
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

We next introduce the variables. u is the x-velocity component while v is the y-velocity component. Note that we use a MONOMIAL basis but we could just as well have used L2_LAGRANGE or L2_HIERARCHIC. Note that the default basis order is 1 or FIRST. Similarly the default family is LAGRANGE so the pressure variable is implicitly LAGRANGE. We note that a MONOMIAL basis does not have the cross-terms that L2_LAGRANGE or L2_HIERARCHIC bases do for non-simplicial elements; consequently for non-simplicial elements it is less expensive to solve while having a worse error constant. Additionally, a MONOMIAL basis has worse conditioning than L2_LAGRANGE or L2_HIERARCHIC although that difference may not be appreciable until high polynomial orders. As a reference in this test case, the condition number with a 5x5 mesh is 741 with a MONOMIAL basis for velocity and 265 for L2_HIERARCHIC and L2_LAGRANGE. We note also that there is no limitation on the family pairings for the hybrid CG-DG scheme so long as the pressure variable is continuous and the velocity variable is discontinuous. Finally, error convergence rates in the L2 norm for an equal order basis of degree nn will be n+1n + 1 for velocity and nn for pressure. The same convergence rates will be observed with a basis of degree nn for velocity and n−1n - 1 for pressure: n+1n + 1 in the L2 norm for velocity and nn in the L2 norm for pressure.

[Variables]
  [u]
    family = MONOMIAL
  []
  [v]
    family = MONOMIAL
  []
  [pressure]
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

We then add the Kernels block which adds the parts of the finite element weak form terms that are integrated over element volumes. The first three kernels comprise the x-momentum equation advection, diffusion, and pressure terms and are of type ADConservativeAdvection, MatDiffusion, and PressureGradient respectively. The next three kernels are the same physics but for the y-momentum equation. The final kernel, mass, reuses the ADConservativeAdvection kernel but is applied to the mass continuity equation which is applied to the pressure variable.

[Kernels]
  [momentum_x_convection]
    type = ADConservativeAdvection
    variable = u
    velocity = 'velocity'
    advected_quantity = 'rhou'
  []
  [momentum_x_diffusion]
    type = MatDiffusion
    variable = u
    diffusivity = 'mu'
  []
  [momentum_x_pressure]
    type = PressureGradient
    integrate_p_by_parts = false
    variable = u
    pressure = pressure
    component = 0
  []
  [momentum_y_convection]
    type = ADConservativeAdvection
    variable = v
    velocity = 'velocity'
    advected_quantity = 'rhov'
  []
  [momentum_y_diffusion]
    type = MatDiffusion
    variable = v
    diffusivity = 'mu'
  []
  [momentum_y_pressure]
    type = PressureGradient
    integrate_p_by_parts = false
    variable = v
    pressure = pressure
    component = 1
  []
  [mass]
    type = ADConservativeAdvection
    variable = pressure
    velocity = velocity
    advected_quantity = -1
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

The DGKernels section adds the weak form terms corresponding to integrations over internal faces. These only exist for the momentum equation since the pressure is discretized with a continuous basis. The first two DGKernels are for the x-momentum equation's advection and diffusion fluxes and are of type ADDGAdvection and DGDiffusion respectively. The last two DGKernels in the block are the same physics added for the y-momentum equation.

[DGKernels]
  [momentum_x_convection]
    type = ADDGAdvection
    variable = u
    velocity = 'velocity'
    advected_quantity = 'rhou'
  []
  [momentum_x_diffusion]
    type = DGDiffusion
    variable = u
    sigma = 6
    epsilon = -1
    diff = 'mu'
  []
  [momentum_y_convection]
    type = ADDGAdvection
    variable = v
    velocity = 'velocity'
    advected_quantity = 'rhov'
  []
  [momentum_y_diffusion]
    type = DGDiffusion
    variable = v
    sigma = 6
    epsilon = -1
    diff = 'mu'
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

The BCs section adds the weak form terms which are integrated over external boundary faces. Diffusive fluxes are applied using the DGFunctionDiffusionDirichletBC object. The desired value of the solution on the boundary is set using the function parameter. In u_walls we weakly impose the condition that the x-velocity (u) is 0 on all but the top boundary. In vu_walls we weakly impose the condition that the y-velocity (v) is 0 on all boundaries. The impact of the lid is applied in u_top where we impose a velocity of ${U} which is 1. Here we point out a final advantage of the DG scheme for velocity is that we can naturally impose discontinuous boundary conditions, e.g. at the corners where the lid meets the walls there is a discontinuity in what we weakly set for the velocity. The final boundary condition in the BCs block, pressure_pin, imposes the constraint the pressure be 0 at the pinned_node.

[BCs]
  [u_walls]
    type = DGFunctionDiffusionDirichletBC
    boundary = 'left bottom right'
    variable = u
    sigma = 6
    epsilon = -1
    function = '0'
    diff = 'mu'
  []
  [v_walls]
    type = DGFunctionDiffusionDirichletBC
    boundary = 'left bottom right top'
    variable = v
    sigma = 6
    epsilon = -1
    function = '0'
    diff = 'mu'
  []
  [u_top]
    type = DGFunctionDiffusionDirichletBC
    boundary = 'top'
    variable = u
    sigma = 6
    epsilon = -1
    function = '${U}'
    diff = 'mu'
  []
  [pressure_pin]
    type = DirichletBC
    variable = pressure
    boundary = 'pinned_node'
    value = 0
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

The Materials block defines properties necessary to complete the simulation. These include the density, rho, the viscosity, mu, the vector velocity, velocity, the x-momentum, rhou, and the y-momentum, rhov.

[Materials]
  [const]
    type = ADGenericConstantMaterial
    prop_names = 'rho'
    prop_values = '${rho}'
  []
  [const_reg]
    type = GenericConstantMaterial
    prop_names = 'mu'
    prop_values = '${mu}'
  []
  [vel]
    type = ADVectorFromComponentVariablesMaterial
    vector_prop_name = 'velocity'
    u = u
    v = v
  []
  [rhou]
    type = ADParsedMaterial
    property_name = 'rhou'
    coupled_variables = 'u'
    material_property_names = 'rho'
    expression = 'rho*u'
  []
  [rhov]
    type = ADParsedMaterial
    property_name = 'rhov'
    coupled_variables = 'v'
    material_property_names = 'rho'
    expression = 'rho*v'
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

The AuxVariables and AuxKernels blocks are added for ease of comparison with MOOSE finite volume results in Paraview.

[AuxVariables]
  [vel_x]
    family = MONOMIAL
    order = CONSTANT
  []
  [vel_y]
    family = MONOMIAL
    order = CONSTANT
  []
  [p]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [vel_x]
    type = ProjectionAux
    variable = vel_x
    v = u
  []
  [vel_y]
    type = ProjectionAux
    variable = vel_y
    v = v
  []
  [p]
    type = ProjectionAux
    variable = p
    v = pressure
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

The simulation execution is defined in the Executioner block. We are performing a steady-state calculation as indicated by type = Steady. Because the Jacobian is completely accurate we use solve_type = NEWTON instead of solve_type = PJFNK. In this case we are solving with a direct solver via -pc_type lu. More sophisticated preconditioning techniques can be constructed using field splits.

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  nl_rel_tol = 1e-12
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

We request Exodus output in the Outputs block

[Outputs]
  exodus = true
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

Finally, we define a ParsedPostprocessor which prints the Reynolds number to console output

[Postprocessors]
  [Re]
    type = ParsedPostprocessor
    expression = '${rho} * ${U} * ${l} / ${mu}'
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_element/ins/cg-dg-hybrid/lid-driven/hybrid-cg-dg.i)

References

  1. Aditya K Pandare and Hong Luo. A hybrid reconstructed discontinuous galerkin and continuous galerkin finite element method for incompressible flows on unstructured grids. Journal of Computational Physics, 322:491–510, 2016.[BibTeX]