circular_driven_cavity.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // Driver for adaptive 2D quarter circle driven cavity. Solved with black
31 // box adaptation, using Taylor Hood and Crouzeix Raviart elements.
32 
33 // Generic oomph-lib header
34 #include "generic.h"
35 
36 // Navier Stokes headers
37 #include "navier_stokes.h"
38 
39 // The mesh
40 #include "meshes/quarter_circle_sector_mesh.h"
41 
42 using namespace std;
43 
44 using namespace oomph;
45 
46 
47 //==start_of_namespace===================================================
48 /// Namespace for physical parameters
49 //=======================================================================
50 namespace Global_Physical_Variables
51 {
52  /// Reynolds number
53  double Re=100;
54 
55  /// Reynolds/Froude number
56  double Re_invFr=100;
57 
58  /// Gravity vector
59  Vector<double> Gravity(2);
60 
61  /// Functional body force
62  void body_force(const double& time, const Vector<double>& x,
63  Vector<double>& result)
64  {
65  result[0]=0.0;
66  result[1]=-Re_invFr;
67  }
68 
69  /// Zero functional body force
70  void zero_body_force(const double& time, const Vector<double>& x,
71  Vector<double>& result)
72  {
73  result[0]=0.0;
74  result[1]=0.0;
75  }
76 
77 } // end_of_namespace
78 
79 //////////////////////////////////////////////////////////////////////
80 //////////////////////////////////////////////////////////////////////
81 //////////////////////////////////////////////////////////////////////
82 
83 //==start_of_problem_class============================================
84 /// Driven cavity problem in quarter circle domain, templated
85 /// by element type.
86 //====================================================================
87 template<class ELEMENT>
88 class QuarterCircleDrivenCavityProblem : public Problem
89 {
90 
91 public:
92 
93  /// Constructor
95  NavierStokesEquations<2>::NavierStokesBodyForceFctPt body_force_fct_pt);
96 
97  /// Destructor: Empty
99 
100  /// Update the after solve (empty)
102 
103  /// \short Update the problem specs before solve.
104  /// (Re-)set velocity boundary conditions just to be on the safe side...
106  {
107  // Setup tangential flow along boundary 0:
108  unsigned ibound=0;
109  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
110  for (unsigned inod=0;inod<num_nod;inod++)
111  {
112  // Tangential flow
113  unsigned i=0;
114  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,1.0);
115  // No penetration
116  i=1;
117  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
118  }
119 
120  // Overwrite with no flow along all other boundaries
121  unsigned num_bound = mesh_pt()->nboundary();
122  for(unsigned ibound=1;ibound<num_bound;ibound++)
123  {
124  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
125  for (unsigned inod=0;inod<num_nod;inod++)
126  {
127  for (unsigned i=0;i<2;i++)
128  {
129  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
130  }
131  }
132  }
133  } // end_of_actions_before_newton_solve
134 
135 
136  /// After adaptation: Unpin pressure and pin redudant pressure dofs.
138  {
139  // Unpin all pressure dofs
140  RefineableNavierStokesEquations<2>::
141  unpin_all_pressure_dofs(mesh_pt()->element_pt());
142 
143  // Pin redundant pressure dofs
144  RefineableNavierStokesEquations<2>::
145  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
146 
147  // Now pin the first pressure dof in the first element and set it to 0.0
148  fix_pressure(0,0,0.0);
149  } // end_of_actions_after_adapt
150 
151  /// Doc the solution
152  void doc_solution(DocInfo& doc_info);
153 
154 private:
155 
156  /// Pointer to body force function
157  NavierStokesEquations<2>::NavierStokesBodyForceFctPt Body_force_fct_pt;
158 
159  /// Fix pressure in element e at pressure dof pdof and set to pvalue
160  void fix_pressure(const unsigned &e, const unsigned &pdof,
161  const double &pvalue)
162  {
163  //Cast to proper element and fix pressure
164  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
165  fix_pressure(pdof,pvalue);
166  } // end_of_fix_pressure
167 
168 }; // end_of_problem_class
169 
170 
171 
172 //==start_of_constructor==================================================
173 /// Constructor for driven cavity problem in quarter circle domain
174 //========================================================================
175 template<class ELEMENT>
177  NavierStokesEquations<2>::NavierStokesBodyForceFctPt body_force_fct_pt) :
178  Body_force_fct_pt(body_force_fct_pt)
179 {
180 
181  // Build geometric object that parametrises the curved boundary
182  // of the domain
183 
184  // Half axes for ellipse
185  double a_ellipse=1.0;
186  double b_ellipse=1.0;
187 
188  // Setup elliptical ring
189  GeomObject* Wall_pt=new Ellipse(a_ellipse,b_ellipse);
190 
191  // End points for wall
192  double xi_lo=0.0;
193  double xi_hi=2.0*atan(1.0);
194 
195  //Now create the mesh
196  double fract_mid=0.5;
197  Problem::mesh_pt() = new
198  RefineableQuarterCircleSectorMesh<ELEMENT>(
199  Wall_pt,xi_lo,fract_mid,xi_hi);
200 
201  // Set error estimator
202  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
203  dynamic_cast<RefineableQuarterCircleSectorMesh<ELEMENT>*>(
204  mesh_pt())->spatial_error_estimator_pt()=error_estimator_pt;
205 
206  // Set the boundary conditions for this problem: All nodes are
207  // free by default -- just pin the ones that have Dirichlet conditions
208  // here: All boundaries are Dirichlet boundaries.
209  unsigned num_bound = mesh_pt()->nboundary();
210  for(unsigned ibound=0;ibound<num_bound;ibound++)
211  {
212  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
213  for (unsigned inod=0;inod<num_nod;inod++)
214  {
215  // Loop over values (u and v velocities)
216  for (unsigned i=0;i<2;i++)
217  {
218  mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
219  }
220  }
221  } // end loop over boundaries
222 
223  //Find number of elements in mesh
224  unsigned n_element = mesh_pt()->nelement();
225 
226  // Loop over the elements to set up element-specific
227  // things that cannot be handled by constructor: Pass pointer to Reynolds
228  // number
229  for(unsigned e=0;e<n_element;e++)
230  {
231  // Upcast from GeneralisedElement to the present element
232  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
233 
234  //Set the Reynolds number, etc
235  el_pt->re_pt() = &Global_Physical_Variables::Re;
236  //Set the Re/Fr
237  el_pt->re_invfr_pt() = &Global_Physical_Variables::Re_invFr;
238  //Set Gravity vector
239  el_pt->g_pt() = &Global_Physical_Variables::Gravity;
240  //set body force function
241  el_pt->body_force_fct_pt() = Body_force_fct_pt;
242 
243  } // end loop over elements
244 
245  // Initial refinement level
246  refine_uniformly();
247  refine_uniformly();
248 
249  // Pin redudant pressure dofs
250  RefineableNavierStokesEquations<2>::
251  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
252 
253  // Now pin the first pressure dof in the first element and set it to 0.0
254  fix_pressure(0,0,0.0);
255 
256  // Setup equation numbering scheme
257  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
258 
259 } // end_of_constructor
260 
261 
262 
263 //==start_of_doc_solution=================================================
264 /// Doc the solution
265 //========================================================================
266 template<class ELEMENT>
268 {
269 
270  ofstream some_file;
271  char filename[100];
272 
273  // Number of plot points
274  unsigned npts=5;
275 
276 
277  // Output solution
278  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
279  doc_info.number());
280  some_file.open(filename);
281  mesh_pt()->output(some_file,npts);
282  some_file.close();
283 
284 } // end_of_doc_solution
285 
286 
287 
288 
289 //==start_of_main======================================================
290 /// Driver for QuarterCircleDrivenCavityProblem test problem
291 //=====================================================================
292 int main()
293 {
294 
295  // Set output directory and initialise count
296  DocInfo doc_info;
297  doc_info.set_directory("RESLT");
298 
299  // Set max. number of black-box adaptation
300  unsigned max_adapt=3;
301 
302  // Solve problem 1 with Taylor-Hood elements
303  //--------------------------------------------
304  {
305  // Set up downwards-Gravity vector
308 
309  // Set up Gamma vector for stress-divergence form
310  NavierStokesEquations<2>::Gamma[0]=1;
311  NavierStokesEquations<2>::Gamma[1]=1;
312 
313  // Build problem with Gravity vector in stress divergence form,
314  // using zero body force function
317 
318  // Solve the problem with automatic adaptation
319  problem.newton_solve(max_adapt);
320 
321  // Step number
322  doc_info.number()=0;
323 
324  // Output solution
325  problem.doc_solution(doc_info);
326 
327  } // end of problem 1
328 
329 
330 
331  // Solve problem 2 with Taylor Hood elements
332  //--------------------------------------------
333  {
334  // Set up zero-Gravity vector
337 
338  // Set up Gamma vector for simplified form
339  NavierStokesEquations<2>::Gamma[0]=0;
340  NavierStokesEquations<2>::Gamma[1]=0;
341 
342  // Build problem with body force function and simplified form,
343  // using body force function
346 
347  // Solve the problem with automatic adaptation
348  problem.newton_solve(max_adapt);
349 
350  // Step number
351  doc_info.number()=1;
352 
353  // Output solution
354  problem.doc_solution(doc_info);
355 
356  } // end of problem 2
357 
358 } // end_of_main
359 
360 
QuarterCircleDrivenCavityProblem(NavierStokesEquations< 2 >::NavierStokesBodyForceFctPt body_force_fct_pt)
Constructor.
Vector< double > Gravity(2)
Gravity vector.
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Functional body force.
void zero_body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Zero functional body force.
void actions_before_newton_solve()
Update the problem specs before solve. (Re-)set velocity boundary conditions just to be on the safe s...
~QuarterCircleDrivenCavityProblem()
Destructor: Empty.
NavierStokesEquations< 2 >::NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
void fix_pressure(const unsigned &e, const unsigned &pdof, const double &pvalue)
Fix pressure in element e at pressure dof pdof and set to pvalue.
double Re
Reynolds number.
void doc_solution(DocInfo &doc_info)
Doc the solution.
int main()
Driver for QuarterCircleDrivenCavityProblem test problem.
double Re_invFr
Reynolds/Froude number.
void actions_after_newton_solve()
Update the after solve (empty)
void actions_after_adapt()
After adaptation: Unpin pressure and pin redudant pressure dofs.