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<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  [gen]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = 0
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = ${l}
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = 0
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = ${l}
    nx<<<{"description": "Number of elements in the X direction"}>>> = 20
    ny<<<{"description": "Number of elements in the Y direction"}>>> = 20
  []
  [corner_node]
    type = ExtraNodesetGenerator<<<{"description": "Creates a new node set and a new boundary made with the nodes the user provides.", "href": "../../source/meshgenerators/ExtraNodesetGenerator.html"}>>>
    new_boundary<<<{"description": "The names of the boundaries to create"}>>> = 'pinned_node'
    nodes<<<{"description": "The nodes you want to be in the nodeset (Either this parameter or \"coord\" must be supplied)."}>>> = '0'
    input<<<{"description": "The mesh we want to modify"}>>> = gen
  []
[]
(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 will be for velocity and for pressure. The same convergence rates will be observed with a basis of degree for velocity and for pressure: in the L2 norm for velocity and in the L2 norm for pressure.

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [u]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
  [v]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
  [pressure]
  []
[]
(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 ConservativeAdvection, MatDiffusion, and PressureGradient respectively. The next three kernels are the same physics but for the y-momentum equation. The final kernel, mass, reuses the ConservativeAdvection kernel but is applied to the mass continuity equation which is applied to the pressure variable.

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [momentum_x_convection]
    type = ADConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../source/kernels/ConservativeAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    velocity = 'velocity'
    advected_quantity<<<{"description": "An optional material property to be advected. If not supplied, then the variable will be used."}>>> = 'rhou'
  []
  [momentum_x_diffusion]
    type = MatDiffusion<<<{"description": "Diffusion equation Kernel that takes an isotropic Diffusivity from a material property", "href": "../../source/kernels/MatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = 'mu'
  []
  [momentum_x_pressure]
    type = PressureGradient<<<{"description": "Implements the pressure gradient term for one of the Navier Stokes momentum equations.", "href": "../../source/kernels/PressureGradient.html"}>>>
    integrate_p_by_parts<<<{"description": "Whether to integrate the pressure term by parts"}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    pressure<<<{"description": "pressure"}>>> = pressure
    component<<<{"description": "number of component (0 = x, 1 = y, 2 = z)"}>>> = 0
  []
  [momentum_y_convection]
    type = ADConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../source/kernels/ConservativeAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = v
    velocity = 'velocity'
    advected_quantity<<<{"description": "An optional material property to be advected. If not supplied, then the variable will be used."}>>> = 'rhov'
  []
  [momentum_y_diffusion]
    type = MatDiffusion<<<{"description": "Diffusion equation Kernel that takes an isotropic Diffusivity from a material property", "href": "../../source/kernels/MatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = v
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = 'mu'
  []
  [momentum_y_pressure]
    type = PressureGradient<<<{"description": "Implements the pressure gradient term for one of the Navier Stokes momentum equations.", "href": "../../source/kernels/PressureGradient.html"}>>>
    integrate_p_by_parts<<<{"description": "Whether to integrate the pressure term by parts"}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = v
    pressure<<<{"description": "pressure"}>>> = pressure
    component<<<{"description": "number of component (0 = x, 1 = y, 2 = z)"}>>> = 1
  []
  [mass]
    type = ADConservativeAdvection<<<{"description": "Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.", "href": "../../source/kernels/ConservativeAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    velocity = velocity
    advected_quantity<<<{"description": "An optional material property to be advected. If not supplied, then the variable will be used."}>>> = -1
  []
[]
(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<<<{"href": "../../syntax/DGKernels/index.html"}>>>]
  [momentum_x_convection]
    type = ADDGAdvection<<<{"description": "Adds internal face advection flux contributions for discontinuous Galerkin discretizations", "href": "../../source/dgkernels/ADDGAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    velocity<<<{"description": "Velocity vector"}>>> = 'velocity'
    advected_quantity<<<{"description": "An optional material property to be advected. If not supplied, then the variable will be used."}>>> = 'rhou'
  []
  [momentum_x_diffusion]
    type = DGDiffusion<<<{"description": "Computes residual contribution for the diffusion operator using discontinous Galerkin method.", "href": "../../source/dgkernels/DGDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    sigma<<<{"description": "sigma"}>>> = 6
    epsilon<<<{"description": "epsilon"}>>> = -1
    diff<<<{"description": "The diffusion (or thermal conductivity or viscosity) coefficient."}>>> = 'mu'
  []
  [momentum_y_convection]
    type = ADDGAdvection<<<{"description": "Adds internal face advection flux contributions for discontinuous Galerkin discretizations", "href": "../../source/dgkernels/ADDGAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = v
    velocity<<<{"description": "Velocity vector"}>>> = 'velocity'
    advected_quantity<<<{"description": "An optional material property to be advected. If not supplied, then the variable will be used."}>>> = 'rhov'
  []
  [momentum_y_diffusion]
    type = DGDiffusion<<<{"description": "Computes residual contribution for the diffusion operator using discontinous Galerkin method.", "href": "../../source/dgkernels/DGDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = v
    sigma<<<{"description": "sigma"}>>> = 6
    epsilon<<<{"description": "epsilon"}>>> = -1
    diff<<<{"description": "The diffusion (or thermal conductivity or viscosity) coefficient."}>>> = 'mu'
  []
[]
(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<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [u_walls]
    type = DGFunctionDiffusionDirichletBC<<<{"description": "Diffusion Dirichlet boundary condition for discontinuous Galerkin method.", "href": "../../source/bcs/DGFunctionDiffusionDirichletBC.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left bottom right'
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    sigma<<<{"description": "Sigma"}>>> = 6
    epsilon<<<{"description": "Epsilon"}>>> = -1
    function<<<{"description": "The forcing function."}>>> = '0'
    diff<<<{"description": "The diffusion (or thermal conductivity or viscosity) coefficient."}>>> = 'mu'
  []
  [v_walls]
    type = DGFunctionDiffusionDirichletBC<<<{"description": "Diffusion Dirichlet boundary condition for discontinuous Galerkin method.", "href": "../../source/bcs/DGFunctionDiffusionDirichletBC.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left bottom right top'
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = v
    sigma<<<{"description": "Sigma"}>>> = 6
    epsilon<<<{"description": "Epsilon"}>>> = -1
    function<<<{"description": "The forcing function."}>>> = '0'
    diff<<<{"description": "The diffusion (or thermal conductivity or viscosity) coefficient."}>>> = 'mu'
  []
  [u_top]
    type = DGFunctionDiffusionDirichletBC<<<{"description": "Diffusion Dirichlet boundary condition for discontinuous Galerkin method.", "href": "../../source/bcs/DGFunctionDiffusionDirichletBC.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'top'
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    sigma<<<{"description": "Sigma"}>>> = 6
    epsilon<<<{"description": "Epsilon"}>>> = -1
    function<<<{"description": "The forcing function."}>>> = '${U}'
    diff<<<{"description": "The diffusion (or thermal conductivity or viscosity) coefficient."}>>> = 'mu'
  []
  [pressure_pin]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'pinned_node'
    value<<<{"description": "Value of the BC"}>>> = 0
  []
[]
(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<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [const]
    type = ADGenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../../source/materials/GenericConstantMaterial.html"}>>>
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = 'rho'
    prop_values<<<{"description": "The values associated with the named properties"}>>> = '${rho}'
  []
  [const_reg]
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../../source/materials/GenericConstantMaterial.html"}>>>
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = 'mu'
    prop_values<<<{"description": "The values associated with the named properties"}>>> = '${mu}'
  []
  [vel]
    type = ADVectorFromComponentVariablesMaterial<<<{"description": "Computes a vector material property from coupled variables", "href": "../../source/materials/VectorFromComponentVariablesMaterial.html"}>>>
    vector_prop_name<<<{"description": "The name to give the declared vector material property"}>>> = 'velocity'
    u<<<{"description": "x-component"}>>> = u
    v<<<{"description": "y-component"}>>> = v
  []
  [rhou]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'rhou'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'u'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'rho'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'rho*u'
  []
  [rhov]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'rhov'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'v'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'rho'
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = 'rho*v'
  []
[]
(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<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [vel_x]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [vel_y]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [p]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [vel_x]
    type = ProjectionAux<<<{"description": "Returns the specified variable as an auxiliary variable with a projection of the source variable. If they are the same type, this amounts to a simple copy.", "href": "../../source/auxkernels/ProjectionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = vel_x
    v<<<{"description": "Variable to take the value of."}>>> = u
  []
  [vel_y]
    type = ProjectionAux<<<{"description": "Returns the specified variable as an auxiliary variable with a projection of the source variable. If they are the same type, this amounts to a simple copy.", "href": "../../source/auxkernels/ProjectionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = vel_y
    v<<<{"description": "Variable to take the value of."}>>> = v
  []
  [p]
    type = ProjectionAux<<<{"description": "Returns the specified variable as an auxiliary variable with a projection of the source variable. If they are the same type, this amounts to a simple copy.", "href": "../../source/auxkernels/ProjectionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = p
    v<<<{"description": "Variable to take the value of."}>>> = pressure
  []
[]
(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<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  nl_rel_tol = 1e-12
[]
(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<<<{"href": "../../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[]
(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<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [Re]
    type = ParsedPostprocessor<<<{"description": "Computes a parsed expression with post-processors", "href": "../../source/postprocessors/ParsedPostprocessor.html"}>>>
    expression<<<{"description": "function expression"}>>> = '${rho} * ${U} * ${l} / ${mu}'
  []
[]
(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]