Using CDO schemes for solving the Navier-Stokes system

Table of Contents

• Introduction
• Settings done in cs_user_model()
• Settings done in cs_user_finalize_setup()
  • Boundary conditions
  • Definition of source terms
• Settings done in cs_user_parameters()
  • Set the strategy to solve the Navier-Stokes system when a monolithic coupling is used
• To go further

Introduction

The Navier–Stokes equations can be solved using CDO face-based discretizations. Other space discretizations are not available. Up to now, there is no turbulence modelling when using CDO schemes. This is a work in progress.

The rationale to set up the computation in the case of the Navier-Stokes equations is the very near of the one explained for user-defined equations here. One only focuses on the specifities arising from the Navier-Stokes case.

Settings done in cs_user_model()

The activation of the NavSto module is done thanks to the function cs_navsto_system_activate For instance, here are two examples to activate and set the main parameters for the NavSto module

  {
    cs_domain_t  *domain = cs_glob_domain;
 
    cs_navsto_system_activate(domain->boundaries,
                              /* Model: type + additional flag */
                              CS_NAVSTO_MODEL_STOKES,
                              CS_NAVSTO_MODEL_STEADY,
                              /* Velocity-pressure coupling */
                              CS_NAVSTO_COUPLING_MONOLITHIC,
                              /* Post flag */
                              CS_NAVSTO_POST_VELOCITY_DIVERGENCE   |
                              0);
  }

or

  {
    cs_domain_t  *domain = cs_glob_domain;
    cs_flag_t  post_flag =
      CS_NAVSTO_POST_VELOCITY_DIVERGENCE |
      CS_NAVSTO_POST_VELOCITY_GRADIENT   |
      CS_NAVSTO_POST_MASS_DENSITY        |
      CS_NAVSTO_POST_PRESSURE_GRADIENT;
 
    cs_navsto_system_activate(domain->boundaries,
                              /* Model: type + additional flag */
                              CS_NAVSTO_MODEL_INCOMPRESSIBLE_NAVIER_STOKES,
                              0,
                              /* Velocity-pressure coupling */
                              CS_NAVSTO_COUPLING_MONOLITHIC,
                              /* Post flag */
                              post_flag);
  }

The first parameter is the structure managing the domain boundaries. An example of settings for the domain boundaries is described here.

The second parameter specifies the type of model to consider among the following choice:

• CS_NAVSTO_MODEL_STOKES,
• CS_NAVSTO_MODEL_OSEEN,
• CS_NAVSTO_MODEL_INCOMPRESSIBLE_NAVIER_STOKES,

The base model can be updated thanks to the third parameter which is a flag built from the following elemental bit (cs_navsto_param_model_bit_t):

• CS_NAVSTO_MODEL_STEADY to specify the model is steady-state (by default, this is not the case).
• CS_NAVSTO_MODEL_GRAVITY_EFFECTS
• CS_NAVSTO_MODEL_BOUSSINESQ

The fourth parameter specifies which type of velocity-pressure algorithm will be used (cs_navsto_param_coupling_t). The choice is done among:

• CS_NAVSTO_COUPLING_ARTIFICIAL_COMPRESSIBILITY (cf. [14] and [15]),
• CS_NAVSTO_COUPLING_MONOLITHIC
• CS_NAVSTO_COUPLING_PROJECTION (work in progress),

The last parameter specifies predefined post-processing operations. As the third parameter, this a flag built from the following elemental bit (cs_navsto_param_post_bit_t):

• CS_NAVSTO_POST_VELOCITY_DIVERGENCE
• CS_NAVSTO_POST_KINETIC_ENERGY
• CS_NAVSTO_POST_VORTICITY
• CS_NAVSTO_POST_VELOCITY_GRADIENT
• CS_NAVSTO_POST_STREAM_FUNCTION (This adds an equation named streamfunction_eq and its variable field named stream_function. This equation relies on a CDO vertex-based equation.)
• CS_NAVSTO_POST_HELICITY
• CS_NAVSTO_POST_ENSTROPHY
• CS_NAVSTO_POST_MASS_DENSITY
• CS_NAVSTO_POST_CELL_MASS_FLUX_BALANCE
• CS_NAVSTO_POST_PRESSURE_GRADIENT

Predefined equations associated to the Navier–Stokes equations are

• the momentum equation is automatically added when activated with a monolithic or artificial compressibility velocity-pressure coupling

In all cases, a vector-valued field named velocity and a scalar-valued field named _"pressure"_ are created. Moreover, several properties are added:

• the property mass_density (a cs_property_t structure);
• the property laminar viscosity (a cs_property_t structure);
• the properties turbulent_viscosity and the total_viscosity are added if the model of turbulence is different from the laminar one (cf. cs_turb_model_t);
• the property graddiv_coef (a cs_property_t structure) when the artificial compressibility is set;

along with the advection field mass_flux (a cs_adv_field_t structure)

Settings done in cs_user_finalize_setup()

Boundary conditions

In the case of the NavSto module, this is done as follows (One does not access to the equation directly since it depends on the way the velocity/pressure coupling is done).

  {
    cs_navsto_param_t  *nsp = cs_navsto_system_get_param();
 
    /* Define the boundary conditions from a function called "_vel_def" */
 
    cs_navsto_set_velocity_inlet_by_analytic(nsp,
                                             "boundary_faces",
                                             _vel_def,
                                             NULL);
  }

where the function _vel_def is defined as follows

static inline void
__vel(const cs_real_3_t pxyz,
      cs_real_3_t       res)
{
  const cs_real_t x  = pxyz[0], y = pxyz[1], z = pxyz[2];
  const cs_real_t sx = sin(_2pi*x), cx = cos(_2pi*x),
                  sy = sin(_2pi*y), cy = cos(_2pi*y),
                  sz = sin(_2pi*z), cz = cos(_2pi*z);
 
  res[0] = .5 * sx * cy * cz;
  res[1] = .5 * cx * sy * cz;
  res[2] = -    cx * cy * sz;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
static void
_vel_def(cs_real_t           time,
         cs_lnum_t           n_pts,
         const cs_lnum_t    *pt_ids,
         const cs_real_t    *xyz,
         bool                dense_output,
         void               *input,
         cs_real_t          *res)
{
  CS_UNUSED(input);
  CS_UNUSED(time);
 
  if (pt_ids != NULL && !dense_output) {
 
#   pragma omp parallel for if(n_pts > CS_THR_MIN)
    for (cs_lnum_t p = 0; p < n_pts; p++) {
      const cs_lnum_t id = pt_ids[p];
      __vel(xyz + 3*id, res + 3*id);
    }
 
  }
  else if (pt_ids != NULL && dense_output) {
 
#   pragma omp parallel for if(n_pts > CS_THR_MIN)
    for (cs_lnum_t p = 0; p < n_pts; p++) {
      const cs_lnum_t id = pt_ids[p];
      __vel(xyz + 3*id, res + 3*p);
    }
 
  }
  else {
 
#   pragma omp parallel for if(n_pts > CS_THR_MIN)
    for (cs_lnum_t p = 0; p < n_pts; p++)
      __vel(xyz + 3*p, res + 3*p);
 
  }
}

Definition of source terms

Since the equation on which this source term applies, depends on the choice of the velocity-pressure algorithm, the way to proceed varies slightly of the way used on a user-defined equation.

  {
    cs_navsto_param_t  *nsp = cs_navsto_system_get_param();
 
    /* Define the source term */
 
    cs_xdef_t *st_def =
      cs_navsto_add_source_term_by_analytic(nsp, "cells", _src_def, NULL);
 
    /* Optionally, the parameters of this definition can be updated */
 
    cs_xdef_set_quadrature(st_def, CS_QUADRATURE_BARY_SUBDIV);
 
  }

where the function _src_def is defined as follows

static inline void
__st(const cs_real_3_t pxyz,
     cs_real_3_t       res)
{
  const cs_real_t x  = pxyz[0], y = pxyz[1], z = pxyz[2];
  const cs_real_t sx = sin(_2pi*x), cx = cos(_2pi*x),
                  sy = sin(_2pi*y), cy = cos(_2pi*y),
                  sz = sin(_2pi*z), cz = cos(_2pi*z);
  res[0] =   6.*_pis * sx * cy * cz + _2pi * cx * sy * sz;
  res[1] =   6.*_pis * cx * sy * cz + _2pi * sx * cy * sz;
  res[2] = -12.*_pis * cx * cy * sz + _2pi * sx * sy * cz;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
static void
_src_def(cs_real_t           time,
         cs_lnum_t           n_pts,
         const cs_lnum_t    *pt_ids,
         const cs_real_t    *xyz,
         bool                dense_output,
         void               *input,
         cs_real_t          *res)
{
  CS_UNUSED(input);
  CS_UNUSED(time);
 
  if (pt_ids != NULL && !dense_output) {
 
#   pragma omp parallel for if(n_pts > CS_THR_MIN)
    for (cs_lnum_t p = 0; p < n_pts; p++) {
      const cs_lnum_t id = pt_ids[p];
      __st(xyz + 3*id, res + 3*id);
    }
 
  }
  else if (pt_ids != NULL && dense_output) {
 
#   pragma omp parallel for if(n_pts > CS_THR_MIN)
    for (cs_lnum_t p = 0; p < n_pts; p++) {
      const cs_lnum_t id = pt_ids[p];
      __st(xyz + 3*id, res + 3*p);
    }
 
  }
  else {
 
#   pragma omp parallel for if(n_pts > CS_THR_MIN)
    for (cs_lnum_t p = 0; p < n_pts; p++)
      __st(xyz + 3*p, res + 3*p);
 
  }
}

Settings done in cs_user_parameters()

The rationale is similar to the one detailed in this section.

Here are listed some specificities related to the NavSto module.

Set the strategy to solve the Navier-Stokes system when a monolithic coupling is used

When a monolithic velocity-pressure is set, the linear system to solve is a saddle-point problem. This class of linear systems needs specific choices of preconditioner/solver. The default settings is not always the optimal choice in terms of efficiency. The settings of saddle-point problems is detailed here

To go further

The detailed description of CDO schemes, the mathmatical analysis and numerical results on benchmarks are available in the following PhD thesis:

• PhD thesis: "Compatible Discrete Operator schemes for the unsteady incompessible Navier-Stokes equations" [15]

Additional publications :

• Article: "Artificial compressibility methods for the incompressible Navier-Stokes equations using lowest-order face-based schemes on polytopal meshes" [14]