Example problem: Adaptive solution of the 2D driven cavity problem

In a previous example we demonstrated the solution of the 2D driven cavity problem using oomph-lib's 2D Taylor-Hood and Crouzeix-Raviart Navier-Stokes elements on a uniform mesh. The computed solution was clearly under-resolved near the corners of the domain where the discontinuity in the velocity boundary conditions creates pressure singularities.

In this example we shall re-solve the driven cavity problem with the refineable versions of oomph-lib's quadrilateral Navier-Stokes elements – the RefineableQTaylorHoodElement<2> and the RefineableQCrouzeixRaviartElement<2>. Enabling spatial adaptivity for this problem involves the same straightforward steps as for a scalar problem:

Two additional steps tend to be required during the adaptive solution of Navier-Stokes problems:

The driver code discussed below illustrates the use of these functions. The section Comments and Exercises provides a more detailed discussion of the technical details involved in their implementation.



The example problem

We shall illustrate the spatially adaptive solution of the steady 2D Navier-Stokes equations by re-considering the 2D steady driven cavity problem:

The 2D steady driven cavity problem in a square domain.
Solve

\[ Re\phantom{i}u_j\frac{\partial u_i}{\partial x_j} = - \frac{\partial p}{\partial x_i} + \frac{\partial }{\partial x_j} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right), \ \ \ \ \ \ \ \ \ \ (1) \]

and

\[ \frac{\partial u_i}{\partial x_i} = 0, \]

in the square domain $ D = \left\{x_i \in [0,1]; i=1,2 \right\} $, subject to the Dirichlet boundary conditions

\[ \left. \mathbf{u}\right|_{\partial D}=(0,0), \ \ \ \ \ \ \ \ \ \ (2) \]

on right, top and left boundaries and

\[ \left. \mathbf{u}\right|_{\partial D}=(1,0), \ \ \ \ \ \ \ \ \ \ (3) \]

on the bottom boundary, $ x_2 = 0 $ .

Solution with Crouzeix-Raviart elements

The figure below shows "carpet plots" of the velocity and pressure fields as well as a contour plot of the pressure distribution with superimposed streamlines for a Reynolds number of $ Re=100$. The velocity vanishes along the entire domain boundary, apart from the bottom boundary $ (x_2 = 0) $ where the moving "lid" imposes a unit tangential velocity which drives a large vortex, centred at $ (x_1,x_2) \approx (0.62,0.26)$. The pressure singularities created by the velocity discontinuities at $ (x_1,x_2)=(0,0) $ and $ (x_1,x_2)=(1,0) $ are now much better resolved.

adaptive_driven_cavity_CR.gif
Plot of the velocity and pressure fields computed with adaptive Crouzeix-Raviart elements for Re=100.

Solution with Taylor-Hood elements

The next figure shows the corresponding results obtained from a computation with adaptive Taylor-Hood elements.

adaptive_driven_cavity_TH.gif
Plot of the velocity and pressure fields computed with adaptive Taylor-Hood elements for Re=100.


Global parameters and functions

The global namespace used to define the problem parameters is identical to the one in the non-adaptive version.

//==start_of_namespace===================================================
/// Namespace for physical parameters
//=======================================================================
namespace Global_Physical_Variables
{
/// Reynolds number
double Re=100;
} // end_of_namespace



The driver code

The main driver code is virtually identical to that in the non-adaptive version. We specify the appropriate refineable element types and use the black-box adaptive Newton solver, allowing for up to three levels of spatial adaptivity.

//==start_of_main======================================================
/// Driver for RefineableDrivenCavity test problem
//=====================================================================
int main()
{
// Set output directory
DocInfo doc_info;
doc_info.set_directory("RESLT");
// Set max. number of black-box adaptation
unsigned max_adapt=3;
// Solve problem with Taylor Hood elements
//---------------------------------------
{
//Build problem
// Solve the problem with automatic adaptation
problem.newton_solve(max_adapt);
// Step number
doc_info.number()=0;
//Output solution
problem.doc_solution(doc_info);
} // end of Taylor Hood elements
// Solve problem with Crouzeix Raviart elements
//--------------------------------------------
{
// Build problem
// Solve the problem with automatic adaptation
problem.newton_solve(max_adapt);
// Step number
doc_info.number()=1;
//Output solution
problem.doc_solution(doc_info);
} // end of Crouzeix Raviart elements
} // end_of_main



The problem class

Most of the problem class is identical to that in the non-adaptive version of the code : We provide a constructor and destructor and use the function Problem::actions_before_newton_solve() to (re-)assign the boundary conditions.

//==start_of_problem_class============================================
/// Driven cavity problem in rectangular domain, templated
/// by element type.
//====================================================================
template<class ELEMENT>
class RefineableDrivenCavityProblem : public Problem
{
public:
/// Constructor
/// Destructor: Empty
/// Update the after solve (empty)
/// \short Update the problem specs before solve.
/// (Re-)set velocity boundary conditions just to be on the safe side...
{
// Setup tangential flow along boundary 0:
unsigned ibound=0;
unsigned num_nod= mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
{
// Tangential flow
unsigned i=0;
mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,1.0);
// No penetration
i=1;
mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
}
// Overwrite with no flow along all other boundaries
unsigned num_bound = mesh_pt()->nboundary();
for(unsigned ibound=1;ibound<num_bound;ibound++)
{
unsigned num_nod= mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
{
for (unsigned i=0;i<2;i++)
{
mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
}
}
}
} // end_of_actions_before_newton_solve

As discussed in the introduction, we use the function Problem::actions_after_adapt() to ensure that, regardless of the mesh adaptation pattern, exactly one pressure degree of freedom is pinned. We start by unpinning all pressure degrees of freedom:

/// After adaptation: Unpin pressure and pin redudant pressure dofs.
void actions_after_adapt()
{
// Unpin all pressure dofs
RefineableNavierStokesEquations<2>::
unpin_all_pressure_dofs(mesh_pt()->element_pt());

[Note that this function (implemented as a static member function of the NavierStokesEquations<DIM> class) unpins the pressure degrees of freedom in all elements that are specified by the input argument (a vector of pointers to the these elements). This implementation allows certain elements in a mesh to be excluded from the procedure; this is required in problems where a mesh contains multiple element types. In the present problem, the mesh contains only Navier-Stokes elements, so we pass a vector of pointers to all elements in the mesh (returned by the function Mesh::element_pt()] to the function.]

Following the mesh adaptation any redundant nodal pressures must be pinned, so that hanging pressure degrees of freedom are treated correctly. We note that calling this function is essential for Taylor-Hood elements. The function may be executed without any adverse effect for all other Navier-Stokes elements; see Comments and Exercises for more details on the implementation.

// Pin redundant pressure dofs
RefineableNavierStokesEquations<2>::
pin_redundant_nodal_pressures(mesh_pt()->element_pt());

Finally, we pin a single pressure degree of freedom (the first pressure value in the first element in the mesh) and set its value to zero.

// Now set the first pressure dof in the first element to 0.0
fix_pressure(0,0,0.0);
} // end_of_actions_after_adapt

The remainder of the problem class remains as before.

/// Doc the solution
void doc_solution(DocInfo& doc_info);
private:
///Fix pressure in element e at pressure dof pdof and set to pvalue
void fix_pressure(const unsigned &e, const unsigned &pdof,
const double &pvalue)
{
//Cast to proper element and fix pressure
dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
fix_pressure(pdof,pvalue);
} // end_of_fix_pressure
}; // end_of_problem_class



The problem constructor

The constructor remains largely as before. We create an adaptive mesh, build and assign an error estimator and pin the redundant nodal pressure degrees of freedom.

//==start_of_constructor==================================================
/// Constructor for RefineableDrivenCavity problem
///
//========================================================================
template<class ELEMENT>
{
// Setup mesh
// # of elements in x-direction
unsigned n_x=10;
// # of elements in y-direction
unsigned n_y=10;
// Domain length in x-direction
double l_x=1.0;
// Domain length in y-direction
double l_y=1.0;
// Build and assign mesh
Problem::mesh_pt() =
new RefineableRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
// Set error estimator
Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
dynamic_cast<RefineableRectangularQuadMesh<ELEMENT>*>(mesh_pt())->
spatial_error_estimator_pt()=error_estimator_pt;
// Set the boundary conditions for this problem: All nodes are
// free by default -- just pin the ones that have Dirichlet conditions
// here: All boundaries are Dirichlet boundaries.
unsigned num_bound = mesh_pt()->nboundary();
for(unsigned ibound=0;ibound<num_bound;ibound++)
{
unsigned num_nod= mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
{
// Loop over values (u and v velocities)
for (unsigned i=0;i<2;i++)
{
mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
}
}
} // end loop over boundaries
//Find number of elements in mesh
unsigned n_element = mesh_pt()->nelement();
// Loop over the elements to set up element-specific
// things that cannot be handled by constructor: Pass pointer to Reynolds
// number
for(unsigned e=0;e<n_element;e++)
{
// Upcast from GeneralisedElement to the present element
ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
//Set the Reynolds number, etc
el_pt->re_pt() = &Global_Physical_Variables::Re;
} // end loop over elements
// Pin redudant pressure dofs
RefineableNavierStokesEquations<2>::
pin_redundant_nodal_pressures(mesh_pt()->element_pt());
// Now set the first pressure dof in the first element to 0.0
fix_pressure(0,0,0.0);
// Setup equation numbering scheme
cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
} // end_of_constructor



Post-processing

The post-processing function is identical to that in the non-adaptive version of the code.

//==start_of_doc_solution=================================================
/// Doc the solution
//========================================================================
template<class ELEMENT>
{
ofstream some_file;
char filename[100];
// Number of plot points
unsigned npts=5;
// Output solution
sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
mesh_pt()->output(some_file,npts);
some_file.close();
} // end_of_doc_solution



Comments and Exercises

Hanging nodes in problems with vector-valued unknowns

We discussed in an earlier example for a scalar (Poisson) problem how oomph-lib's mesh adaptation routines create hanging nodes and how the values that are stored at such nodes are automatically constrained to ensure the inter-element continuity of the solution. The methodology employed for scalar problems is easily generalised to problems with vector-valued unknowns, provided that all unknowns are represented by the same isoparametric interpolation between the elements' nodal values. In such problems, the unknown nodal values at the hanging nodes are constrained to be linear combination of the corresponding values at their "master nodes". The list of "master nodes" and the corresponding "hanging weights" (contained in a HangInfo object) are the same for all unknowns.

To allow the use spatial adaptivity for elements in which different unknowns are represented by different interpolation schemes (e.g. in 2D quadrilateral Taylor-Hood elements where the two velocity components are represented by bi-quadratic interpolation between the values stored at the element's 3x3 nodes, whereas the pressure is represented by bi-linear interpolation between the pressure values stored at the element's 2x2 corner nodes) oomph-lib allows the different nodal values to have their own list of "master nodes" and "hanging weights". This is achieved as follows;

Adaptivity for Taylor-Hood elements

The issues

In non-adaptive 2D [3D] Taylor-Hood elements, every node stores (at least) two [three] nodal values which represent the two [three] velocity components. The four corner [eight vertex] nodes store an additional third [fourth] value which represents the pressure. If the mesh is subjected to non-uniform refinement, some of the mid-side nodes in large elements also act as corner nodes for adjacent smaller elements, as illustrated in the figure below.

pressure_hang.gif
An adapted mesh, with nodes storing a pressure value denoted by circles, and velocity nodes by squares.

The figure illustrates that the "hanging status" of the various degrees of freedom can be become fairly involved. For instance

To illustrate that oomph-lib's automatic mesh adaptation procedures are able to deal with these cases, the figure below shows a "carpet plot" of the pressure distribution, $ p(x_1,x_2) $, obtained from a (strongly under-resolved) driven-cavity computation on the mesh shown above. The figure illustrates that the hanging node constraints ensure the inter-element continuity of the pressure throughout the domain, and, in particular, along the "right" boundary of the largest element.

hanging_pressure.gif
The under-resolved pressure distribution for the driven cavity problem, using the above mesh and computed with Taylor-Hood elements.

Details of the implementation

To facilitate the book-keeping for such problems, all nodes in the refineable Taylor-Hood elements store three [four] nodal values, even though, depending on the mesh's refinement pattern some of the pressure values will not be used. To eliminate the "redundant" pressure degrees of freedom from the problem, we provide the function

NavierStokesEquations<DIM>::pin_redundant_nodal_pressures(...)

which pins the "redundant" pressure degrees of freedom in all elements specified by the input argument (a vector of pointers to the Navier-Stokes elements). This function must be called after the initial mesh has been created, and after each mesh adaptation. The function first pins all nodal pressure values, using the function

NavierStokesEquations<DIM>::pin_all_nodal_pressure_dofs()

and then unpins the nodal pressure values at the elements' corner [vertex] nodes, using the function

NavierStokesEquations<DIM>::unpin_proper_nodal_pressure_dofs()

These functions are implemented as empty virtual functions in the NavierStokesEquations<DIM> class which provides a base class for all Navier-Stokes elements. The empty functions are overwritten for QTaylorHoodElement<DIM> and remain empty for all other Navier-Stokes elements, therefore they can be called for any element type.

Adaptivity for Crouzeix-Raviart elements

As discussed in the previous example, oomph-lib's isoparametric 2D [3D] Crouzeix-Raviart elements employ a piecewise bi- [tri-]linear, globally discontinuous pressure representation. In each element, the pressure is represented by bi-[tri-]linear basis functions, multiplied by 3 [4] pressure values which are stored in the element's internal Data. Since the pressure representation is discontinuous, the pressure values do not have to be subjected to any constraints to ensure inter-element continuity. Each Node stores two [three] velocity degrees of freedom. Since the velocity representation is isoparametric, the default assignment for the nodal values' HangInfo pointer is appropriate and no further action is required.

Exercises

  1. Confirm that a warning message is issued if the function NavierStokesEquations<DIM>::pin_all_nodal_pressure_dofs() is not called following the mesh adaptation.
  2. Investigate how the pressure distribution changes with each adaptation. [Hint: You can call doc_solution(...) from actions_after_newton_solve() to document the progress of the mesh adaptation.]


Source files for this tutorial



PDF file

A pdf version of this document is available.