No Matches
COR-10: Navier-Stokes equation


Navier-Stokes equations (NSE), governing the motion of a viscous fluid, are used in various applications: from simulations of the flow in blood vessels to studies of the air flow around aeroplane wings and rotor blades, scaling up to models of ocean and atmospheric currents. Even in the case of an incompressible steady flow, NSE are non-linear due to the effect of the inertia, which is more pronounced in case of a higher Reynolds number. In this tutorial we discuss the implementation of a viscous fluid model in MoFEM using hierarchical basis functions. This approach permits us to locally increase the order of approximation, enforcing conformity across finite element boundaries, without the need to change the implementation of an element. Moreover, the requirement of different approximation orders for primal (velocity) and dual (pressure) variables, necessary for a simulation of the flow using the mixed formulation, can be easily satisfied.

Problem statement

An incompressible isoviscous steady-state flow in a domain \(\Omega\) is governed by the following equations:

\[ \begin{align} \label{eq:balance_momentum} \rho \left(\mathbf{u}\cdot\nabla\right)\mathbf{u} - \mu\nabla^2\mathbf{u} + \nabla p &= \mathbf{f},\\ \label{eq:cont} \nabla\cdot\mathbf{u} &= 0, \end{align} \]

where \eqref{eq:balance_momentum} is the set of Navier-Stokes equations, representing the balance of the momentum, and \eqref{eq:cont} is the continuity equation; \(\mathbf{u}=[u_x, u_y, u_z]^\intercal\) is the velocity field, \(p\) is the hydrostatic pressure field, \(\rho\) is the fluid mass density, \(\mu\) is fluid viscosity and $\mathbf{f}$ is the density of external forces. The boundary value problem complements the above equations by the Dirichlet and Neumann conditions on the boundary \(\partial\Omega\):

\[ \begin{align} \label{eq:bc_d}\mathbf{u} = \mathbf{u}_D & \;\text{on}\; \Gamma_D,\\ \label{eq:bc_n}\mathbf{n}\cdot \left[-p\mathbf{I}+\mu\left(\nabla\mathbf{u} + \nabla\mathbf{u}^\intercal\right)\right] = \mathbf{g}_N & \;\text{on} \;\Gamma_N, \end{align} \]

where \(\mathbf{u}_D\) is the prescribed velocity on the part of the boundary \(\Gamma_D\subset\partial\Omega\), and \(\mathbf{g}_N\) is the prescribed traction vector on the part of the boundary \(\Gamma_N\subset\partial\Omega\), \(\mathbf{n}\) is an outward normal.

In this tutorial we will solve the problem of the fluid flow around a rigid sphere (see Figure 1) of a radius \(r\) positioned in the centre of a cubic domain, the side length \(2l\) of which is considered sufficiently large compared to \(r\), so that a uniform far-field velocity on the exterior boundaries is valid (note that the body forces are neglected). Exploiting the symmetry of the problem, we will use a quarter of the domain in our simulations.

Figure 1: (a) Finite-element mesh used for simulation of fluid flow around a rigid sphere. (b) Sketch of the problem set-up on a section $z=0$ of the mesh.

Note that due to the structure of differential operators in the Navier-Stokes equations \eqref{eq:balance_momentum}, the fluid pressure has to be specified on a part of boundary in order to obtain a unique solution. To achieve that, we will combine Dirichlet and Neumann condition on the outlet boundary of the domain, see see Figure 1. In the case when the flow is perpendicular to the boundary ( \(u_y = u_z = 0\)), which is enforced by the Dirichlet boundary condition, the hydrostatic fluid pressure equals to the normal traction on such boundary, see Theorem 1 in [10].

Non-dimensionalization and scaling

Reynolds number is introduced for the problems governed by the Navier-Stokes equations as a measure of the ratio of inertial forces to viscous forces:

\[ \begin{equation} \mathcal{R} = \frac{\rho U L }{\mu}, \end{equation} \]

where \(U\) is the scale for the velocity and \(L\) is a relevant length scale. For the problem of the fluid flow around a sphere, considered in this tutorial, the far-field velocity plays the role of the velocity scale, while \(L=2r\), where \(r\) is the radius of the sphere.

When viscous forces dominate over the inertia forces, \(\mathcal{R} \ll 1\), the non-linear term in \eqref{eq:balance_momentum} can be neglected, simplifying NSE down to Stokes equations. However, often this is not the case, i.e. the nonlinear term may become dominant over the viscous forces. In these cases non-dimensionalization (scaling) of Navier-Stokes equations is very useful, which permits to decrease the number of coefficients to only one – Reynolds number \(\mathcal{R}\). The following scales can be used for the case of flow with relatively low Reynolds number (e.g. less than 100).

Non-dimensionalization of Navier-Stokes equations
Physical quantity Scale Dimensionless variable
Length \(L\) \(\mathbf{x}^{*}=\mathbf{x}\:/\:L, \quad \nabla^{*}(\cdot) = L\nabla(\cdot)\)
Velocity \(U\) \(\mathbf{u}^{*}=\mathbf{u}\:/\:U\)
Pressure \(P=\frac{\mu U}{L}\) \(p^{*}=p\:/\:P\)

Mote that the scale for pressure depends on the scales for the length and the velocity. Using these scales, Navier-stokes equations (1) - (2) can be written in the following dimensionless form (considering zero external forces):

\[ \begin{align} \label{eq:balance_momentum_dim_less} \mathcal{R} \left(\mathbf{u}^{*}\cdot\nabla^{*}\right)\mathbf{u}^{*} - (\nabla^{*})^2\mathbf{u}^{*} + \nabla^{*} p^{*} &= 0,\\ \label{eq:cont_dim_less} \nabla^{*}\cdot\mathbf{u}^{*} &= 0, \end{align} \]

Note that now the influence of the non-linearity (inertia terms) on the problem is fully controlled by the value of the Reynolds number. On the one hand, if \(\mathcal{R} \ll 1\), or if one simply wants to consider the Stokes equation, the nonlinear term can be dropped. On the other hand, if the Reynolds number is sufficiently high, then the non-linearity can become strong and may pose problems for convergence (see below discussion of the linearisation of the problem for the Newton-Raphson method). In this case a certain iterative procedure can be used to obtain a solution. Indeed, to find the solution for a given Reynolds number, a number of intermediate problems can be solved for range of smaller \(\mathcal{R}\), eventually reaching the original value. Such technique will also be used in this tutorial. Note also that below we will drop the * symbol, implying that during computation all considered quantities are dimensionless, and a transformation to the dimensional variables is performed before post-processing of the results takes place.

Finite-element formulation

Weak form

The weak statement of the problem \eqref{eq:balance_momentum}-\eqref{eq:bc_n} reads: Find a vector field \(\mathbf{u} \in \mathbf{H}^1(\Omega)\) and a scalar field \(p \in L^2(\Omega)\), such that for any test functions \(\mathbf{v} \in \mathbf{H}^1_0(\Omega)\) and \(q \in L^2(\Omega)\):

\[ \begin{align} \label{eq:weak} \mathcal{R}\int\limits_\Omega \left(\mathbf{u}\cdot\nabla\right)\mathbf{u} \cdot\mathbf{v} \,d\Omega + \int\limits_{\Omega}\nabla\mathbf{u}\mathbin{:}\nabla\mathbf{v}\, d\Omega - \int\limits_{\Omega}p\, \nabla\cdot\mathbf{v} \, d\Omega &= \int\limits_{\Omega}\mathbf{f}\cdot\mathbf{v}\,d\Omega + \int\limits_{\Gamma_N}\mathbf{g}_N\cdot\mathbf{v}\,d\Gamma_N,\\ - \int\limits_{\Omega}q\, \nabla\cdot\mathbf{u} \, d\Omega &= 0 \end{align} \]

Upon finite-element discretization of the domain \(\Omega\), we consider interpolation of both unknown fields introducing shape functions on each element: \begin{equation} \label{eq:shape} u_i = \sum\limits_{\alpha=1}^{n_{\mathbf{u}}} N_{\alpha}\, u_i^\alpha, \quad p = \sum\limits_{\beta=1}^{n_p} \Phi_{\beta}\, p^\beta; \quad v_i = \sum\limits_{\alpha=1}^{n_{\mathbf{u}}} N_{\alpha}\, v_i^\alpha, \quad q = \sum\limits_{\beta=1}^{n_p} \Phi_{\beta}\, q^\beta, \end{equation} where \(n_{\mathbf{u}}\) is the number of shape functions associated with the velocity field, and $n_p$ is the similar number for the pressure field. Using the hierarchical basis approximation, the vector of the shape functions can be decomposed into four sub-vectors, consisting of shape functions associated with element's entities: vertices, edges, faces and the volume of the element, e.g. for the velocity field:

\[ \begin{equation} \label{eq:shape_hier} \mathbf{N}^\textit{el} = \left[N_1, \ldots N_\alpha, \ldots N_{n_{\mathbf{u}}}\right]^\intercal = \left[\mathbf{N}^\textit{ver}, \mathbf{N}^\textit{edge}, \mathbf{N}^\textit{face}, \mathbf{N}^\textit{vol}\right]^\intercal. \end{equation} \]

Note on the choice of spaces for the mixed problem

Note that in the weak statement velocity and pressure fields belong to different spaces: \(H^1\) and \(L^2\), respectively. Two different approaches can be used for implementing such mixed formulation in the finite-element framework. One can use generalised Taylor-Hood elements, which feature continuous ( \(H^1\) ) approximation of both velocity and pressure fields, however, in order to enforce stability, the approximations functions for the pressure field should be one order lower than those for the velocity field. Furthermore, Taylor-Hood elements impose certain constraints on the mesh, see [52]. Alternatively, a discontinuous approximation for pressure can be used, however, in this case the difference between the orders of approximation functions for velocity and for pressure should be 2.


The example class has necessary fields to store input parameters and internal data structures, as well as a set of functions for setting up and solving the problem:

commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
CHKERR SNESDestroy(&snes);
CHKERR DMDestroy(&dM);
MoFEMErrorCode runProblem();
PetscBool isPartitioned = PETSC_FALSE;
int orderVelocity = 2; // default approximation orderVelocity
int orderPressure = 1; // default approximation orderPressure
int numHoLevels = 0;
PetscBool isDiscontPressure = PETSC_FALSE;
PetscBool isHoGeometry = PETSC_TRUE;
double pressureScale = 1.0;
double lengthScale = 1.0;
double velocityScale = 1.0;
double reNumber = 1.0;
double density;
double viscosity;
PetscBool isStokesFlow = PETSC_FALSE;
int numSteps = 1;
double lambdaStep = 1.0;
double lambda = 0.0;
int step;
BitRefLevel bitLevel;
DM dM;
SNES snes;
boost::shared_ptr<NavierStokesElement::CommonData> commonData;
SmartPetscObj<Vec> D;
SmartPetscObj<Vec> F;
SmartPetscObj<Mat> Aij;
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feLhsPtr;
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feRhsPtr;
boost::shared_ptr<FaceElementForcesAndSourcesCore> feDragPtr;
boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feDragSidePtr;
boost::shared_ptr<DirichletDisplacementBc> dirichletBcPtr;
boost::ptr_map<std::string, NeumannForcesSurface> neumannForces;
boost::shared_ptr<PostProcVol> fePostProcPtr;
boost::shared_ptr<PostProcFace> fePostProcDragPtr;
MoFEMErrorCode readInput();
MoFEMErrorCode findBlocks();
MoFEMErrorCode setupFields();
MoFEMErrorCode defineFiniteElements();
MoFEMErrorCode setupDiscreteManager();
MoFEMErrorCode setupAlgebraicStructures();
MoFEMErrorCode setupElementInstances();
MoFEMErrorCode setupSNES();
MoFEMErrorCode solveProblem();
MoFEMErrorCode postProcess();
#define CHKERR
Inline error check.
Definition: definitions.h:535
Deprecated interface functions.
[Example Navier Stokes]
MoFEMErrorCode setupElementInstances()
[Setup algebraic structures]
MoFEMErrorCode setupDiscreteManager()
[Define finite elements]
boost::shared_ptr< FaceElementForcesAndSourcesCore > feDragPtr
MoFEMErrorCode defineFiniteElements()
[Setup fields]
boost::shared_ptr< PostProcFace > fePostProcDragPtr
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhsPtr
MoFEMErrorCode setupFields()
[Find blocks]
MoFEMErrorCode setupSNES()
[Setup element instances]
MoFEMErrorCode runProblem()
[Example Navier Stokes]
SmartPetscObj< Vec > F
boost::shared_ptr< NavierStokesElement::CommonData > commonData
boost::shared_ptr< DirichletDisplacementBc > dirichletBcPtr
PetscBool isDiscontPressure
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > feDragSidePtr
MoFEMErrorCode findBlocks()
[Read input]
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhsPtr
boost::ptr_map< std::string, NeumannForcesSurface > neumannForces
MoFEM::Interface & mField
boost::shared_ptr< PostProcVol > fePostProcPtr
MoFEMErrorCode solveProblem()
[Setup SNES]
MoFEMErrorCode postProcess()
[Solve problem]
MoFEMErrorCode readInput()
[Run problem]
SmartPetscObj< Vec > D
SmartPetscObj< Mat > Aij

The set of functions used in this tutorials in the following:

#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416

The function ExampleNavierStokes::runProblem is executed from the main function:

int main(int argc, char *argv[]) {
const char param_file[] = "param_file.petsc";
// Initialise MoFEM
try {
// Create mesh database
moab::Core mb_instance; // create database
moab::Interface &moab = mb_instance; // create interface to database
// Create MoFEM database and link it to MoAB
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
NavierStokesExample ex(m_field);
CHKERR ex.runProblem();
// finish work cleaning memory, getting statistics, etc
return 0;
std::string param_file
static char help[]
int main()
Definition: adol-c_atom.cpp:46
Catch errors.
Definition: definitions.h:372
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112

Reading mesh and input parameters

The workflow of solving a problem using the finite-element method starts from reading the mesh and other input parameters:

char mesh_file_name[255];
PetscBool flg_mesh_file;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "NAVIER_STOKES problem",
CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_mesh_file);
CHKERR PetscOptionsBool("-my_is_partitioned", "is partitioned?", "",
isPartitioned, &isPartitioned, PETSC_NULL);
CHKERR PetscOptionsInt("-my_order_u", "approximation orderVelocity", "",
orderVelocity, &orderVelocity, PETSC_NULL);
CHKERR PetscOptionsInt("-my_order_p", "approximation orderPressure", "",
orderPressure, &orderPressure, PETSC_NULL);
CHKERR PetscOptionsInt("-my_num_ho_levels",
"number of higher order boundary levels", "",
numHoLevels, &numHoLevels, PETSC_NULL);
CHKERR PetscOptionsBool("-my_discont_pressure", "discontinuous pressure", "",
CHKERR PetscOptionsBool("-my_ho_geometry", "use second order geometry", "",
isHoGeometry, &isHoGeometry, PETSC_NULL);
CHKERR PetscOptionsScalar("-my_length_scale", "length scale", "", lengthScale,
&lengthScale, PETSC_NULL);
CHKERR PetscOptionsScalar("-my_velocity_scale", "velocity scale", "",
velocityScale, &velocityScale, PETSC_NULL);
CHKERR PetscOptionsBool("-my_stokes_flow", "stokes flow", "", isStokesFlow,
&isStokesFlow, PETSC_NULL);
CHKERR PetscOptionsInt("-my_step_num", "number of steps", "", numSteps,
&numSteps, PETSC_NULL);
ierr = PetscOptionsEnd();
if (flg_mesh_file != PETSC_TRUE) {
// Read whole mesh or part of it if partitioned
if (isPartitioned == PETSC_TRUE) {
// This is a case of distributed mesh and algebra. In that case each
// processor keeps only part of the problem.
const char *option;
CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
} else {
// In this case, we have distributed algebra, i.e. assembly of vectors and
// matrices is in parallel, but whole mesh is stored on all processors.
// snes and matrix scales well, however problem set-up of problem is
// not fully parallel.
const char *option;
option = "";
CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 3,
static PetscErrorCode ierr
char mesh_file_name[255]
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

Once the mesh is read, we find the set of tetrahedra corresponding to the computational domain. Furthermore, we identify the fluid-structure interface, i.e. the set of triangles corresponding to the surface of the sphere:

if (bit->getName().compare(0, 9, "MAT_FLUID") == 0) {
const int id = bit->getMeshsetId();
CHKERR mField.get_moab().get_entities_by_dimension(
bit->getMeshset(), 3, commonData->setOfBlocksData[id].eNts, true);
std::vector<double> attributes;
if (attributes.size() < 2) {
"should be at least 2 attributes but is %d",
commonData->setOfBlocksData[id].iD = id;
commonData->setOfBlocksData[id].fluidViscosity = attributes[0];
commonData->setOfBlocksData[id].fluidDensity = attributes[1];
viscosity = commonData->setOfBlocksData[id].fluidViscosity;
density = commonData->setOfBlocksData[id].fluidDensity;
if (bit->getName().compare(0, 9, "INT_SOLID") == 0) {
Range tets, tet;
const int id = bit->getMeshsetId();
CHKERR mField.get_moab().get_entities_by_type(
bit->getMeshset(), MBTRI, commonData->setOfFacesData[id].eNts, true);
CHKERR mField.get_moab().get_adjacencies(
commonData->setOfFacesData[id].eNts, 3, true, tets,
tet = Range(tets.front(), tets.front());
for (auto &bit : commonData->setOfBlocksData) {
if (bit.second.eNts.contains(tet)) {
commonData->setOfFacesData[id].fluidViscosity =
commonData->setOfFacesData[id].fluidDensity = bit.second.fluidDensity;
commonData->setOfFacesData[id].iD = id;
if (commonData->setOfFacesData[id].fluidViscosity < 0) {
"Cannot find a fluid block adjacent to a given solid face");
Definition: definitions.h:148
Definition: definitions.h:31
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
auto bit
set bit

Setting up fields, finite elements and the problem itself

Once the domain and the fluid-structure interface are identified, we can setup the problem. First, we define the velocity and pressure fields. Note that for velocity we use space H1, while for pressure the user can choose between using also H1 space (continuous approximation, Taylor-Hood element), or a discontinuous approximation using L2 space. Note that in the former case, the approximation functions for pressure has to be one order lower than that of the velocity (e.g. order 2 for velocity and order 1 for pressure), while the the latter case the difference between the two has to be 2 orders (e.g. order 3 for velocity and order 1 for pressure).

Furthermore, the properties of hierarchical basis functions can be used to locally increase the approximation order for the elements on the fluid-structure interface, where the gradients of the velocity can be expected to be the highest due to the no-slip boundary condition on the surface of the sphere.

} else {
CHKERR mField.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
CHKERR mField.set_field_order(0, MBVERTEX, "VELOCITY", 1);
CHKERR mField.set_field_order(0, MBEDGE, "VELOCITY", orderVelocity);
CHKERR mField.set_field_order(0, MBTRI, "VELOCITY", orderVelocity);
CHKERR mField.set_field_order(0, MBTET, "VELOCITY", orderVelocity);
CHKERR mField.set_field_order(0, MBVERTEX, "PRESSURE", 1);
CHKERR mField.set_field_order(0, MBEDGE, "PRESSURE", orderPressure);
CHKERR mField.set_field_order(0, MBTRI, "PRESSURE", orderPressure);
CHKERR mField.set_field_order(0, MBTET, "PRESSURE", orderPressure);
if (numHoLevels > 0) {
std::vector<Range> levels(numHoLevels);
Range ents;
for (int ll = 0; ll != numHoLevels; ll++) {
Range verts;
CHKERR mField.get_moab().get_connectivity(ents, verts, true);
for (auto d : {1, 2, 3}) {
CHKERR mField.get_moab().get_adjacencies(verts, d, false, ents,
levels[ll] = subtract(ents, ents.subset_by_type(MBVERTEX));
for (int ll = numHoLevels - 1; ll >= 1; ll--) {
levels[ll] = subtract(levels[ll], levels[ll - 1]);
int add_order = 1;
for (int ll = numHoLevels - 1; ll >= 0; ll--) {
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
CHKERR mField.set_field_order(levels[ll], "VELOCITY",
orderVelocity + add_order);
CHKERR mField.set_field_order(levels[ll], "PRESSURE",
orderPressure + add_order);
CHKERR mField.set_field_order(levels[ll].subset_by_type(MBTET),
"PRESSURE", orderPressure + add_order);
CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
// Set 2nd order of approximation of geometry
if (isHoGeometry) {
Range ents, edges;
CHKERR mField.get_moab().get_entities_by_dimension(0, 3, ents);
CHKERR mField.get_moab().get_adjacencies(ents, 1, false, edges,
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(edges);
CHKERR mField.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
Projection10NodeCoordsOnField ent_method_material(mField,
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.

After that the finite elements operating on the set above fields can be defined. In particular, we define elements for solving Navier-Stokes equations, computation of the drag traction and force on the surface of the surface, and, finally, for computing the contribution of Neumann (natural) boundary conditions.

// Add finite element (this defines element, declaration comes later)
CHKERR NavierStokesElement::addElement(mField, "NAVIER_STOKES", "VELOCITY",
CHKERR NavierStokesElement::addElement(mField, "DRAG", "VELOCITY", "PRESSURE",
"MESH_NODE_POSITIONS", 2, &solidFaces);
// setup elements for loading
// build finite elements
// build adjacencies between elements and degrees of freedom
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const string element_name, const string velocity_field_name, const string pressure_field_name, const string mesh_field_name, const int dim=3, Range *ents=nullptr)
Setting up elements.

Once the finite elements are defined, we will need to create element instances and push necessary operators to their pipelines. First, we setup an element for solving Navier-Stokes (or if one prefers, Stokes) equations. Note that right-hand side and left-hand side elements are considered separately. In case of Navier-Stokes equations, operators are pushed to elements pipelines using function NavierStokesElement::setNavierStokesOperators. In particular, the following operators are pushed to right-hand side:

while for the left-hand side we push these operators:

Upon that, we setup in the same way elements for computing the drag traction and the drag force acting on the sphere using operators NavierStokesElement::OpCalcDragTraction and NavierStokesElement::OpCalcDragForce, respectively. Additionally, operators for the Neumann boundary conditions are set using the function MetaNeumannForces::setMomentumFluxOperators and the element for the Dirichlet boundary conditions is created as well. Finally, operators for post-processing output are created: in the volume and on the fluid-solid interface.

feLhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
feRhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feLhsPtr, true, false, false,
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feRhsPtr, true, false, false,
feDragPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
new VolumeElementForcesAndSourcesCoreOnSide(mField));
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feDragSidePtr, true, false, false,
if (isStokesFlow) {
feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
} else {
feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
"PRESSURE", commonData);
dirichletBcPtr = boost::shared_ptr<DirichletDisplacementBc>(
new DirichletDisplacementBc(mField, "VELOCITY", Aij, D, F));
dirichletBcPtr->methodsOp.push_back(new NavierStokesElement::LoadScale());
// dirichletBcPtr->snes_ctx = FEMethod::CTX_SNESNONE;
dirichletBcPtr->snes_x = D;
// Assemble pressure and traction forces
// for postprocessing:
fePostProcPtr = boost::make_shared<PostProcVol>(mField);
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fePostProcPtr, true, false, false,
auto v_ptr = boost::make_shared<MatrixDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto pos_ptr = boost::make_shared<MatrixDouble>();
auto p_ptr = boost::make_shared<VectorDouble>();
new OpCalculateVectorFieldValues<3>("VELOCITY", v_ptr));
new OpCalculateVectorFieldGradient<3, 3>("VELOCITY", grad_ptr));
new OpCalculateScalarFieldValues("PRESSURE", p_ptr));
new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
using OpPPMap = OpPostProcMapInMoab<3, 3>;
new OpPPMap(
fePostProcPtr->getPostProcMesh(), fePostProcPtr->getMapGaussPts(),
{{"PRESSURE", p_ptr}},
{{"VELOCITY", v_ptr}, {"MESH_NODE_POSITIONS", pos_ptr}},
{{"VELOCITY_GRAD", grad_ptr}},
fePostProcDragPtr = boost::make_shared<PostProcFace>(mField);
new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
new OpPPMap(
{{"MESH_NODE_POSITIONS", pos_ptr}},
fePostProcDragPtr, feDragSidePtr, "NAVIER_STOKES", "VELOCITY", "PRESSURE",
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Post post-proc data at points from hash maps.
Set integration rule to volume elements.
static MoFEMErrorCode setCalcDragOperators(boost::shared_ptr< FaceElementForcesAndSourcesCore > dragFe, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for calculating drag force on the solid surface.
static MoFEMErrorCode setPostProcDragOperators(boost::shared_ptr< PostProcFace > postProcDragPtr, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for post processing output of drag traction.
static MoFEMErrorCode setNavierStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Navier-Stokes equations.
static MoFEMErrorCode setStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Stokes equations.

Note that the class for enforcing Dirichlet boundary conditions DirichletDisplacementBc used in the snippet above requires the algebraic data structures for storing components of the right-hand side vector, left-hand side matrix and the solution vector. Therefore, these data structures need to be created beforehand using PETSc:

D = smartCreateDMVector(dM);
F = smartVectorDuplicate(D);
Aij = smartCreateDMMatrix(dM);
CHKERR VecZeroEntries(F);
CHKERR VecZeroEntries(D);
CHKERR MatZeroEntries(Aij);

However, one may note that creation of PETSCc vectors and matrices requires the instance of the MoFEM Discrete Manager (DM), which needs a certain setup as well:

DMType dm_name = "DM_NAVIER_STOKES";
// Register DM problem
CHKERR DMRegister_MoFEM(dm_name);
CHKERR DMSetType(dM, dm_name);
// Create DM instance
CHKERR DMMoFEMCreateMoFEM(dM, &mField, dm_name, bitLevel);
// Configure DM form line command options s
CHKERR DMSetFromOptions(dM);
// Add elements to dM (only one here)
if (mField.check_finite_element("PRESSURE_FE"))
CHKERR DMMoFEMSetIsPartitioned(dM, isPartitioned);
// setup the DM
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.

Finally, the last preparatory step required is the setup of the SNES functionality:

boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
for (; mit != neumannForces.end(); mit++) {
CHKERR DMMoFEMSNESSetFunction(dM, mit->first.c_str(),
&mit->second->getLoopFe(), NULL, NULL);
boost::shared_ptr<FEMethod> null_fe;
CHKERR DMMoFEMSNESSetFunction(dM, "NAVIER_STOKES", feRhsPtr, null_fe,
CHKERR DMMoFEMSNESSetFunction(dM, DM_NO_ELEMENT, null_fe, null_fe,
CHKERR DMMoFEMSNESSetJacobian(dM, DM_NO_ELEMENT, null_fe, dirichletBcPtr,
CHKERR DMMoFEMSNESSetJacobian(dM, "NAVIER_STOKES", feLhsPtr, null_fe,
CHKERR DMMoFEMSNESSetJacobian(dM, DM_NO_ELEMENT, null_fe, null_fe,
SnesCtx *snes_ctx;
// create snes nonlinear solver
CHKERR DMMoFEMGetSnesCtx(dM, &snes_ctx);
CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
CHKERR SNESSetFromOptions(snes);
Definition: DMMoFEM.hpp:10

Solving the problem and post-processing the results

auto scale_problem = [&](double U, double L, double P) {
CHKERR mField.getInterface<FieldBlas>()->fieldScale(L,
ProjectionFieldOn10NodeTet ent_method_on_10nodeTet(
mField, "MESH_NODE_POSITIONS", true, true);
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
ent_method_on_10nodeTet.setNodes = false;
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
CHKERR mField.getInterface<FieldBlas>()->fieldScale(U, "VELOCITY");
CHKERR mField.getInterface<FieldBlas>()->fieldScale(P, "PRESSURE");
CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
1.0 / pressureScale);
step = 0;
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Viscosity: %6.4e\n", viscosity);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Density: %6.4e\n", density);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Length scale: %6.4e\n", lengthScale);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Velocity scale: %6.4e\n",
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Pressure scale: %6.4e\n",
if (isStokesFlow) {
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: 0 (Stokes flow)\n");
} else {
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: %6.4e\n",
while (lambda < 1.0 - 1e-12) {
if (isStokesFlow) {
reNumber = 0.0;
for (auto &bit : commonData->setOfBlocksData) {
bit.second.inertiaCoef = 0.0;
bit.second.viscousCoef = 1.0;
} else {
for (auto &bit : commonData->setOfBlocksData) {
bit.second.inertiaCoef = reNumber;
bit.second.viscousCoef = 1.0;
CHKERR PetscPrintf(
"Step: %d | Lambda: %6.4e | Inc: %6.4e | Re number: %6.4e \n", step,
CHKERR VecAssemblyBegin(D);
CHKERR VecAssemblyEnd(D);
CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
1.0 / pressureScale);
static Index< 'L', 3 > L
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:523
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:534
string out_file_name;
CHKERR DMoFEMLoopFiniteElements(dM, "NAVIER_STOKES", fePostProcPtr);
std::ostringstream stm;
stm << "out_" << step << ".h5m";
out_file_name = stm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n",
CHKERR fePostProcPtr->writeFile(out_file_name.c_str());
CHKERR VecZeroEntries(commonData->pressureDragForceVec);
CHKERR VecZeroEntries(commonData->shearDragForceVec);
CHKERR VecZeroEntries(commonData->totalDragForceVec);
CHKERR DMoFEMLoopFiniteElements(dM, "DRAG", feDragPtr);
auto get_vec_data = [&](auto vec, std::array<double, 3> &data) {
CHKERR VecAssemblyBegin(vec);
CHKERR VecAssemblyEnd(vec);
const double *array;
CHKERR VecGetArrayRead(vec, &array);
if (mField.get_comm_rank() == 0) {
for (int i : {0, 1, 2})
data[i] = array[i];
CHKERR VecRestoreArrayRead(vec, &array);
std::array<double, 3> pressure_drag;
std::array<double, 3> shear_drag;
std::array<double, 3> total_drag;
CHKERR get_vec_data(commonData->pressureDragForceVec, pressure_drag);
CHKERR get_vec_data(commonData->shearDragForceVec, shear_drag);
CHKERR get_vec_data(commonData->totalDragForceVec, total_drag);
if (mField.get_comm_rank() == 0) {
MOFEM_LOG_C("SELF", Sev::inform,
"Pressure drag force: [ %6.4e, %6.4e, %6.4e ]",
pressure_drag[0], pressure_drag[1], pressure_drag[2]);
MOFEM_LOG_C("SELF", Sev::inform,
"Shear drag force: [ %6.4e, %6.4e, %6.4e ]", shear_drag[0],
shear_drag[1], shear_drag[2]);
MOFEM_LOG_C("SELF", Sev::inform,
"Total drag force: [ %6.4e, %6.4e, %6.4e ]", total_drag[0],
total_drag[1], total_drag[2]);
std::ostringstream stm;
stm << "out_drag_" << step << ".h5m";
out_file_name = stm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "out file %s\n",
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:574
FTensor::Index< 'i', SPACE_DIM > i
char out_file_name[255]
virtual int get_comm_rank() const =0

Running the example

We will consider two cases: first, flow governed by Stokes equation, i.e. neglecting the nonlinear terms, and then we will solve full set Navier-Stokes equations. For both cases we will use the same mesh. We will start by partitioning this mesh:

mofem_part -my_file sphere_test.cub -my_nparts 4 -output_file sphere_test_4.h5m

while mofem_part can be found in $HOME/mofem_install/um/build_release/tools/. Here we partitioned the mesh in 4 parts in order to be executed on 4 cores in a distributed-memory parallel environment. Note that any number of parts less or equal to number of physically available cores can be used.

Stokes equation

mpirun -np 4 ./navier_stokes -my_file sphere_test_4.h5m \
-my_order_u 2 \
-my_order_p 1 \
-my_discont_pressure 0 \
-my_num_ho_levels 1 \
-my_ho_geometry 1 \
-my_step_num 1 \
-my_velocity_scale 0.1 \
-my_length_scale 2 \
-my_stokes_flow 1 \
-my_is_partitioned 1

Navier-Stokes equations

mpirun -np 4 ./navier_stokes -my_file sphere_test_4.h5m \
-my_order_u 2 \
-my_order_p 1 \
-my_discont_pressure 0 \
-my_num_ho_levels 1 \
-my_ho_geometry 1 \
-my_step_num 2 \
-my_velocity_scale 0.1 \
-my_length_scale 2 \
-my_stokes_flow 0 \
-my_is_partitioned 1