Demo problem: Solution of a "free-boundary" Poisson problem -- a simple model for "fluid"-structure interaction.

In this example we shall consider our first (toy!) interaction problem. The problem combines two single-physics problems, studied in earlier examples, and combines them into a coupled free-boundary problem.

We will now consider the coupled problem obtained by using the solution of Poisson's equation at the control node, $ u_{ctrl} $, as the "load", $ f $, that acts on the two circular arcs that define the curvilinear boundaries of $ D_{fish} $. The resulting coupled problem is sketched in the figure below. While this problem is obviously somewhat artificial, it has many of the key features that arise in genuine fluid-structure interaction problems. In particular, the displacement of the domain boundary is driven by the solution of the "bulk" (Poisson) equations, just as the deformation of an elastic structure in an FSI problem is driven by the fluid pressure and shear stresses, i.e. quantities that are derived from the solution of the Navier-Stokes equations in the "bulk" domain.

fish_fsi_sketch.gif
Sketch of the two individual single-physics problems (top) and the coupled problem (bottom).

The two single-physics problems involve two uncoupled sets of equations and unknowns:

The coupling between the two single-physics problem introduces additional dependencies:

We note that most of the methodology required to solve this coupled problem is already available:

The only interaction that still has to be incorporated into the problem formulation is the dependence of the Poisson element's residual vectors on the geometric Data in the ElasticallySupportedRingElement. This interaction arises through the MacroElement/Domain - based node-update function which translates changes in the GeomObject's geometric Data into changes in the nodal positions. Such dependencies may be added to any existing element by "wrapping" the element into the templated wrapper class MacroElementNodeUpdateElement which has the following inheritance structure:

template<class ELEMENT>
class MacroElementNodeUpdateElement : public virtual ELEMENT,
public virtual MacroElementNodeUpdateElementBase

An element of type MacroElementNodeUpdateElement<ELEMENT> is an element of type ELEMENT, and inherits the additional functionality provided by the MacroElementNodeUpdateElementBase base class. The most important additional functionality provided by this class is the ability to add the values stored in the geometric Data of associated GeomObjects to the element's list of unknowns. Once added, the derivatives of the element's residual vector with respect to these additional unknowns are automatically included into the element's Jacobian matrix. This is achieved by overloading the ELEMENT::get_jacobian(...) function and evaluating the additional derivatives by finite differencing. See Comments for details on the implementation.

The solution of the coupled problem therefore only requires a few trivial changes to the single-physics (Poisson) code:

  1. The element type used for the solution of the "bulk" equations must be changed to its "wrapped" counterpart, as discussed above. For instance, if the single-physics code used a nine-node refineable Poisson element of type RefineableQPoissonElement<2,3>, the coupled problem must be discretised by elements of type MacroElementNodeUpdateElement<RefineableQPoissonElement<2,3> > (Yes, it's a bit of a mouthful...).
  2. The "bulk" mesh must be "upgraded" (again via multiple inheritance) to a Mesh that is derived from the MacroElementNodeUpdateMesh base class.
  3. A vector of pointers to those GeomObjects that are involved in an element's MacroElement/Domain - based node update operation must be passed to the elements. (This is done most easily in the constructor of the "upgraded" mesh.) The geometric Data contained in these GeomObjects is then automatically included in the elements' list of unknowns.
  4. The Mesh's node_update() function must be executed whenever the Newton method has changed the values of the unknowns: This is because changing a value that is stored in a GeomObject's geometric Data does not automatically update the positions of any dependent nodes. This is done most easily be including the node_update() function into the Problem::actions_before_newton_convergence_check() function; we refer to another document for a more detailed discussion of the order in which the various "action" functions are called by oomph-lib's Newton solver.

Results

The animation below shows the results of a spatially-adaptive solution of Poisson's equations in the fish-shaped domain, for a variety of domain "heights". This animation was produced with the single-physics Poisson solver discussed in an earlier example.

elastic_fish.gif
Spatially adaptive solution of Poisson's equation in a fish-shaped domain for various `widths' of the domain.

An increase in the height of the domain increases the amplitude of the solution. This is reflected by the red line in the figure below which shows a plot of $ u_{ctrl} $ as a function of $ Y_c $. The green marker shows the solution of the coupled problem for a spring stiffness of $ k=1 $. For this value of the spring stiffness, the solution of the coupled problem should be (and indeed is) located at the intersection of the curve $ u_{ctrl}(Y_c) $ with the diagonal, $ u_{ctrl} = Y_c $, shown by the dashed blue line.

trace.gif
Solution of Poisson's equation at a control node as a function of the `height' of the domain.

Implementation in oomph-lib

The sections below provide the usual annotated listing of the driver code. We stress that only a few trivial changes are required to incorporate the presence of the free boundary into the existing single-physics code:

Global parameters and functions

The namespace ConstSourceForPoisson defines the constant source function, exactly as in the corresponding single-physics code.

//=============start_of_namespace=====================================
/// Namespace for const source term in Poisson equation
//====================================================================
namespace ConstSourceForPoisson
{
/// Const source function
void get_source(const Vector<double>& x, double& source)
{
source = -1.0;
}
} // end of namespace

The Mesh

Meshes that are to be used with MacroElementNodeUpdateElements should be derived (typically by multiple inheritance) from the MacroElementNodeUpdateMesh class. This class overloads the generic Mesh::node_update() function and ensures that the node update is performed by calling the node_update() function of the Mesh's constituent nodes, rather than simply updating their positions, using the FiniteElement::get_x(...) function. The overloaded version is not only more efficient but also ensures that any auxiliary node update functions (e.g. functions that update the no-slip condition on a moving fluid node on a solid boundary) are performed too.

In our driver code we add the additional functionality provided by the MacroElementNodeUpdateMesh class to the RefineableFishMesh class used in the single-physics Poisson problem considered earlier.

//==========start_of_mesh=================================================
/// Refineable, fish-shaped mesh with MacroElement-based node update.
//========================================================================
template<class ELEMENT>
class MyMacroElementNodeUpdateRefineableFishMesh :
public virtual RefineableFishMesh<ELEMENT>,
public virtual MacroElementNodeUpdateMesh
{

The constructor calls the constructors of the underlying RefineableFishMesh. [Note the explicit call to the FishMesh constructor prior to calling the constructor of the RefineableFishMesh. Without this call, only the default (argument-free) constructor of the FishMesh would be called! Consult your favourite C++ book to check on constructors for derived classes if you don't understand this. We recommend Daoqi Yang's brilliant book C++ and Object-Oriented Numeric Computing for Scientists and Engineers.)

public:
/// \short Constructor: Pass pointer to GeomObject that defines
/// the fish's back and pointer to timestepper
/// (defaults to (Steady) default timestepper defined in the Mesh
/// base class).
MyMacroElementNodeUpdateRefineableFishMesh(GeomObject* back_pt,
TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper) :
FishMesh<ELEMENT>(back_pt,time_stepper_pt),
RefineableFishMesh<ELEMENT>(time_stepper_pt)
{

To activate the MacroElementNodeUpdateElement's ability to automatically compute the derivatives of the residual vectors with respect to the geometric Data that determines its nodal positions, we must pass the pointers to the GeomObjects that are involved in the element's MacroElement - based node-update to the elements. In general, an element's node-update will be affected by multiple GeomObjects therefore the set_node_update_info(...) function expects a vector of pointers to GeomObjects. In the present example, only a single GeomObject (the GeomObject that represents the fish's curved "back") determines the nodal position of all elements:

// Set up all the information that's required for MacroElement-based
// node update: Tell the elements that their geometry depends on the
// fishback geometric object.
unsigned n_element = this->nelement();
for(unsigned i=0;i<n_element;i++)
{
// Upcast from FiniteElement to the present element
ELEMENT *el_pt = dynamic_cast<ELEMENT*>(this->element_pt(i));
// There's just one GeomObject
Vector<GeomObject*> geom_object_pt(1);
geom_object_pt[0] = back_pt;
// Tell the element which geom objects its macro-element-based
// node update depends on
el_pt->set_node_update_info(geom_object_pt);
}
} //end of constructor

The destructor can remain empty but we provide a final overload for the Mesh's node_update() function to avoid any ambiguities as to which one is to be used.

/// \short Destructor: empty
virtual ~MyMacroElementNodeUpdateRefineableFishMesh(){}
/// \short Resolve mesh update: Node update current nodal
/// positions via sparse MacroElement-based update.
//void node_update()
// {
// MacroElementNodeUpdateMesh::node_update();
// }
}; // end of mesh class

The driver code

The driver code is very simple: We build the problem with the "wrapped" version of the refineable quadrilateral nine-node Poisson element. Since the initial mesh is very coarse we perform two uniform mesh refinements before solving the problem with automatic spatial adaptivity, allowing for up to two further mesh adaptations.

//==================start_of_main=========================================
/// Driver for "free-boundary" fish poisson solver with adaptation.
//========================================================================
int main()
{
// Shorthand for element type
typedef MacroElementNodeUpdateElement<RefineableQPoissonElement<2,3> >
ELEMENT;
// Build problem
FreeBoundaryPoissonProblem<ELEMENT> problem;
// Do some uniform mesh refinement first
problem.refine_uniformly();
problem.refine_uniformly();
// Solve/doc fully coupled problem, allowing for up to two spatial
// adaptations.
unsigned max_solve=2;
problem.newton_solve(max_solve);
problem.doc_solution();
} // end of main

The problem class

Apart from a few trivial additions, the problem class is virtually identical to that used in the single-physics Poisson problem. The most important addition to the single-physics problem class is the function Problem::actions_before_newton_convergence_check() which updates the nodal positions in the "bulk" Poisson mesh following an update of the geometric Data that controls the position of the curvilinear domain boundary; we refer to another document for a more detailed discussion of the order in which the various "action" functions are called by oomph-lib's Newton solver.

//==========start_of_problem_class====================================
/// Refineable "free-boundary" Poisson problem in deformable
/// fish-shaped domain. Template parameter identifies the element.
//====================================================================
template<class ELEMENT>
class FreeBoundaryPoissonProblem : public Problem
{
public:
/// \short Constructor
FreeBoundaryPoissonProblem();
/// Destructor (empty)
virtual ~FreeBoundaryPoissonProblem(){};
/// Update the problem specs before solve (empty)
void actions_before_newton_solve() {}
/// Update the problem specs after solve (empty)
void actions_after_newton_solve(){}
/// Access function for the fish mesh
MyMacroElementNodeUpdateRefineableFishMesh<ELEMENT>* fish_mesh_pt()
{
return Fish_mesh_pt;
}
/// Doc the solution
void doc_solution();
/// \short Before checking the new residuals in Newton's method
/// we have to update nodal positions in response to possible
/// changes in the position of the domain boundary
void actions_before_newton_convergence_check()
{
fish_mesh_pt()->node_update();
}
private:
/// Pointer to fish mesh
MyMacroElementNodeUpdateRefineableFishMesh<ELEMENT>* Fish_mesh_pt;
/// Pointer to single-element mesh that stores the GeneralisedElement
/// that represents the fish's back
Mesh* Fish_back_mesh_pt;
}; // end of problem class

The Problem constructor

We start by creating the GeomObject/GeneralisedElement that will represent the unknown curvilinear domain boundary and pass it (in its role as a GeomObject) to the constructor of the bulk mesh. We then add the pointer to the bulk mesh to the Problem's collection of submeshes and create an error estimator for the adaptive solution of the Poisson equation.

//=========start_of_constructor===========================================
/// Constructor for adaptive free-boundary Poisson problem in
/// deformable fish-shaped domain.
//========================================================================
template<class ELEMENT>
FreeBoundaryPoissonProblem<ELEMENT>::FreeBoundaryPoissonProblem()
{
// Set coordinates and radius for the circle that will become the fish back
double x_c=0.5;
double y_c=0.0;
double r_back=1.0;
// Build geometric object that will become the fish back
ElasticallySupportedRingElement* fish_back_pt=
new ElasticallySupportedRingElement(x_c,y_c,r_back);
// Build fish mesh with geometric object that specifies the fish back
Fish_mesh_pt=new
MyMacroElementNodeUpdateRefineableFishMesh<ELEMENT>(fish_back_pt);
// Add the fish mesh to the problem's collection of submeshes:
add_sub_mesh(Fish_mesh_pt);
// Create/set error estimator for the fish mesh
fish_mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;

Next we store the pointer to the ElasticallySupportedRingElement in its own Mesh and add it to the Problem's collection of submeshes before building the Problem's global Mesh from its two submeshes:

// Build mesh that will store only the geometric wall element
Fish_back_mesh_pt=new Mesh;
// So far, the mesh is completely empty. Let's add the
// GeneralisedElement that represents the shape
// of the fish's back to it:
Fish_back_mesh_pt->add_element_pt(fish_back_pt);
// Add the fish back mesh to the problem's collection of submeshes:
add_sub_mesh(Fish_back_mesh_pt);
// Now build global mesh from the submeshes
build_global_mesh();

We choose the central node in the Poisson mesh as the control node and use it (in its role as Data) as the "load" for the ElasticallySupportedRingElement.

// Choose a control node: We'll use the
// central node that is shared by all four elements in
// the base mesh because it exists at all refinement levels.
// How many nodes does element 0 have?
unsigned nnod=fish_mesh_pt()->finite_element_pt(0)->nnode();
// The central node is the last node in element 0:
Node* control_node_pt=fish_mesh_pt()->finite_element_pt(0)->node_pt(nnod-1);
// Use the solution (value 0) at the control node as the load
// that acts on the ring. [Note: Node == Data by inheritance]
dynamic_cast<ElasticallySupportedRingElement*>(Fish_mesh_pt->fish_back_pt())->
set_load_pt(control_node_pt);

Finally, we pin the nodal values on all boundaries, apply the homogeneous Dirichlet boundary conditions, pass the pointer to the source function to the elements, and set up the equation numbering scheme.

// Set the boundary conditions for this problem: All nodes are
// free by default -- just pin the ones that have Dirichlet conditions
// here. Set homogeneous boundary conditions everywhere
unsigned num_bound = fish_mesh_pt()->nboundary();
for(unsigned ibound=0;ibound<num_bound;ibound++)
{
unsigned num_nod= fish_mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
{
fish_mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
fish_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,0.0);
}
}
/// Loop over elements and set pointers to source function
unsigned n_element = fish_mesh_pt()->nelement();
for(unsigned i=0;i<n_element;i++)
{
// Upcast from FiniteElement to the present element
ELEMENT *el_pt = dynamic_cast<ELEMENT*>(fish_mesh_pt()->element_pt(i));
//Set the source function pointer
el_pt->source_fct_pt() = &ConstSourceForPoisson::get_source;
}
// Do equation numbering
cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
} // end of constructor

Post-processing

The post-processing routine writes the computed result to an output file.

//============start_of_doc================================================
/// Doc the solution in tecplot format.
//========================================================================
template<class ELEMENT>
void FreeBoundaryPoissonProblem<ELEMENT>::doc_solution()
{
// Number of plot points in each coordinate direction.
unsigned npts=5;
// Output solution
ofstream some_file("RESLT/soln0.dat");
fish_mesh_pt()->output(some_file,npts);
some_file.close();
} // end of doc

Comments

A more detailed description of the theory and the implementation can be found in the paper

and in this talk:

The following subsections provide a brief description of the main features.

Sparse node updates

The key feature of our implementation which allows the efficient computation of the "shape derivatives" is the ability of MacroElementNodeUpdateNodes (discussed in more detail below) to "update their own position" in response to changes in shape/position of the domain boundary. This capability is demonstrated in the following simple example code.

We start by building the Mesh as before

//==================start_of_main=========================================
/// Driver to document sparse MacroElement-based node update.
//========================================================================
int main()
{
// Shorthand for element type
typedef MacroElementNodeUpdateElement<RefineableQPoissonElement<2,3> >
ELEMENT;
// Set coordinates and radius for the circle that will become the fish back
double x_c=0.5;
double y_c=-0.2;
double r_back=1.0;
// Build geometric object that will become the fish back
ElasticallySupportedRingElement* Fish_back_pt=
new ElasticallySupportedRingElement(x_c,y_c,r_back);
// Build fish mesh with geometric object that specifies the fish back
MacroElementNodeUpdateRefineableFishMesh<ELEMENT>* Fish_mesh_pt=new
MacroElementNodeUpdateRefineableFishMesh<ELEMENT>(Fish_back_pt);

and document the mesh (i.e. the shape of its constituent finite elements and the nodal positions):

// Number of plot points in each coordinate direction.
unsigned npts=11;
ofstream some_file;
char filename[100];
// Output initial mesh
unsigned count=0;
sprintf(filename,"RESLT/soln%i.dat",count);
some_file.open(filename);
Fish_mesh_pt->output(some_file,npts);
some_file.close();
count++;

Next, we "manually" increment $ Y_c $, i.e. the y-coordinate of the centre of the circular arc that defines the upper curvilinear boundary of the fish mesh.

// Increment y_c
Fish_back_pt->y_c()+=0.2;

This step mimics the incrementation of one of the Problems's unknowns (recall that in the free-boundary problem considered above, $ Y_c $ has to be determined as part of the solution!) during the finite-difference based computation of the shape derivatives.

For meshes that are not derived from the MacroElementNodeUpdateMeshBase class, the only way to update the nodal positions in response to a change in the boundary position, is to call the Mesh::node_update() function. This updates the position of all nodes in the mesh – a very costly operation.

Meshes that are derived from the MacroElementNodeUpdateMeshBase class contain MacroElementNodeUpdateNodes which can update their own position, as shown here:

// Adjust each node in turn and doc
unsigned nnod=Fish_mesh_pt->nnode();
for (unsigned i=0;i<nnod;i++)
{
// Update individual nodal position
Fish_mesh_pt->node_pt(i)->node_update();
// Doc mesh
sprintf(filename,"RESLT/soln%i.dat",count);
some_file.open(filename);
Fish_mesh_pt->output(some_file,npts);
some_file.close();
count++;
}
} // end of main

We note that the Node::node_update() function is defined as an empty virtual function in the Node base class, indicating that "normal" Nodes cannot "update their own position". The function is overloaded in the MacroElementNodeUpdateNode class, details of which are given below. Overloaded versions of this function also exist in various other derived Node classes (such as as the AlgebraicNodes and the SpineNodes) for which algebraic node update operations are defined.

Here is an animation that illustrates how the successive update of the individual nodal positions in response to the change in the boundary position gradually updates the entire mesh.

sparse_node_update.gif
Illustration of the sparse node-update procedure.

How it works

The implementation employs three key components:

The method also works for non-"toy" problems!

The above example demonstrated how easy it is to "upgrade" a driver code for the solution of a single-physics problem to a fluid-structure-interaction-like free-boundary problem. It is important to stress that the methodology employed in our "toy" free-boundary problem can also be used for genuine fluid-structure interaction problems. For instance, the driver code for the simulation of 2D unsteady finite-Reynolds number flow in a channel with an oscillating wall whose motion is prescribed can easily be extended to a driver code for the corresponding fluid-structure interaction problem in which the wall is replaced by a flexible membrane that is loaded by the fluid traction.

Source files for this tutorial

PDF file

A pdf version of this document is available.