circular_driven_cavity2.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_Variables2
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 QuarterCircleDrivenCavityProblem2 : 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 1:
108  unsigned ibound=1;
109  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
110  for (unsigned inod=0;inod<num_nod;inod++)
111  {
112  // get coordinates
113  double x=mesh_pt()->boundary_node_pt(ibound,inod)->x(0);
114  double y=mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
115  // find Lagrangian coordinate (the angle)
116  double zeta=0.0;
117  if (x!=0.0)
118  {
119  zeta=atan(y/x);
120  }
121  // Tangential flow u0
122  unsigned i=0;
123  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,-sin(zeta));
124  // Tangential flow u1
125  i=1;
126  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,cos(zeta));
127  }
128 
129  // Overwrite with no flow along all boundaries
130  unsigned num_bound = mesh_pt()->nboundary();
131  for(unsigned ibound=0;ibound<num_bound;ibound++)
132  {
133  if (ibound!=1)
134  {
135  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
136  for (unsigned inod=0;inod<num_nod;inod++)
137  {
138  for (unsigned i=0;i<2;i++)
139  {
140  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
141  }
142  }
143  }
144  }
145 
146  } // end_of_actions_before_newton_solve
147 
148 
149  /// After adaptation: Unpin pressure and pin redudant pressure dofs.
151  {
152  // Unpin all pressure dofs
153  RefineableNavierStokesEquations<2>::
154  unpin_all_pressure_dofs(mesh_pt()->element_pt());
155 
156  // Pin redundant pressure dofs
157  RefineableNavierStokesEquations<2>::
158  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
159 
160  // Now pin the first pressure dof in the first element and set it to 0.0
161  fix_pressure(0,0,0.0);
162  } // end_of_actions_after_adapt
163 
164  /// Doc the solution
165  void doc_solution(DocInfo& doc_info);
166 
167 private:
168 
169  /// Pointer to body force function
170  NavierStokesEquations<2>::NavierStokesBodyForceFctPt Body_force_fct_pt;
171 
172  /// Fix pressure in element e at pressure dof pdof and set to pvalue
173  void fix_pressure(const unsigned &e, const unsigned &pdof,
174  const double &pvalue)
175  {
176  //Cast to proper element and fix pressure
177  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
178  fix_pressure(pdof,pvalue);
179  } // end_of_fix_pressure
180 
181 }; // end_of_problem_class
182 
183 
184 
185 //==start_of_constructor==================================================
186 /// Constructor for driven cavity problem in quarter circle domain
187 //========================================================================
188 template<class ELEMENT>
190  NavierStokesEquations<2>::NavierStokesBodyForceFctPt body_force_fct_pt) :
191  Body_force_fct_pt(body_force_fct_pt)
192 {
193 
194  // Build geometric object that parametrises the curved boundary
195  // of the domain
196 
197  // Half axes for ellipse
198  double a_ellipse=1.0;
199  double b_ellipse=1.0;
200 
201  // Setup elliptical ring
202  GeomObject* Wall_pt=new Ellipse(a_ellipse,b_ellipse);
203 
204  // End points for wall
205  double xi_lo=0.0;
206  double xi_hi=2.0*atan(1.0);
207 
208  //Now create the mesh
209  double fract_mid=0.5;
210  Problem::mesh_pt() = new
211  RefineableQuarterCircleSectorMesh<ELEMENT>(
212  Wall_pt,xi_lo,fract_mid,xi_hi);
213 
214  // Set error estimator
215  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
216  dynamic_cast<RefineableQuarterCircleSectorMesh<ELEMENT>*>(
217  mesh_pt())->spatial_error_estimator_pt()=error_estimator_pt;
218 
219  // Set the boundary conditions for this problem: All nodes are
220  // free by default -- just pin the ones that have Dirichlet conditions
221  // here: All boundaries are Dirichlet boundaries.
222  unsigned num_bound = mesh_pt()->nboundary();
223  for(unsigned ibound=0;ibound<num_bound;ibound++)
224  {
225  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
226  for (unsigned inod=0;inod<num_nod;inod++)
227  {
228  // Loop over values (u and v velocities)
229  for (unsigned i=0;i<2;i++)
230  {
231  mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
232  }
233  }
234  } // end loop over boundaries
235 
236  //Find number of elements in mesh
237  unsigned n_element = mesh_pt()->nelement();
238 
239  // Loop over the elements to set up element-specific
240  // things that cannot be handled by constructor: Pass pointer to Reynolds
241  // number
242  for(unsigned e=0;e<n_element;e++)
243  {
244  // Upcast from GeneralisedElement to the present element
245  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
246 
247  //Set the Reynolds number, etc
248  el_pt->re_pt() = &Global_Physical_Variables2::Re;
249  //Set the Re/Fr
250  el_pt->re_invfr_pt() = &Global_Physical_Variables2::Re_invFr;
251  //Set Gravity vector
252  el_pt->g_pt() = &Global_Physical_Variables2::Gravity;
253  //set body force function
254  el_pt->body_force_fct_pt() = Body_force_fct_pt;
255 
256  } // end loop over elements
257 
258  // Initial refinement level
259  refine_uniformly();
260  refine_uniformly();
261 
262  // Pin redudant pressure dofs
263  RefineableNavierStokesEquations<2>::
264  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
265 
266  // Now pin the first pressure dof in the first element and set it to 0.0
267  fix_pressure(0,0,0.0);
268 
269  // Setup equation numbering scheme
270  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
271 
272 } // end_of_constructor
273 
274 
275 
276 //==start_of_doc_solution=================================================
277 /// Doc the solution
278 //========================================================================
279 template<class ELEMENT>
281 {
282 
283  ofstream some_file;
284  char filename[100];
285 
286  // Number of plot points
287  unsigned npts=5;
288 
289 
290  // Output solution
291  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
292  doc_info.number());
293  some_file.open(filename);
294  mesh_pt()->output(some_file,npts);
295  some_file.close();
296 
297 } // end_of_doc_solution
298 
299 
300 
301 
302 //==start_of_main======================================================
303 /// Driver for QuarterCircleDrivenCavityProblem2 test problem
304 //=====================================================================
305 int main()
306 {
307 
308  // Set output directory and initialise count
309  DocInfo doc_info;
310  doc_info.set_directory("RESLT2");
311 
312  // Set max. number of black-box adaptation
313  unsigned max_adapt=3;
314 
315  // Solve problem 1 with Taylor-Hood elements
316  //--------------------------------------------
317  {
318  // Set up downwards-Gravity vector
321 
322  // Set up Gamma vector for stress-divergence form
323  NavierStokesEquations<2>::Gamma[0]=1;
324  NavierStokesEquations<2>::Gamma[1]=1;
325 
326  // Build problem with Gravity vector in stress divergence form,
327  // using zero body force function
330 
331  // Solve the problem with automatic adaptation
332  problem.newton_solve(max_adapt);
333 
334  // Step number
335  doc_info.number()=0;
336 
337  // Output solution
338  problem.doc_solution(doc_info);
339 
340  } // end of problem 1
341 
342 
343 
344  // Solve problem 2 with Taylor Hood elements
345  //--------------------------------------------
346  {
347  // Set up zero-Gravity vector
350 
351  // Set up Gamma vector for simplified form
352  NavierStokesEquations<2>::Gamma[0]=0;
353  NavierStokesEquations<2>::Gamma[1]=0;
354 
355  // Build problem with body force function and simplified form,
356  // using body force function
359 
360  // Solve the problem with automatic adaptation
361  problem.newton_solve(max_adapt);
362 
363  // Step number
364  doc_info.number()=1;
365 
366  // Output solution
367  problem.doc_solution(doc_info);
368 
369  } // end of problem 2
370 
371 } // end_of_main
372 
373 
void actions_after_newton_solve()
Update the after solve (empty)
void body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Functional body force.
QuarterCircleDrivenCavityProblem2(NavierStokesEquations< 2 >::NavierStokesBodyForceFctPt body_force_fct_pt)
Constructor.
int main()
Driver for QuarterCircleDrivenCavityProblem2 test problem.
Vector< double > Gravity(2)
Gravity vector.
NavierStokesEquations< 2 >::NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
void actions_before_newton_solve()
Update the problem specs before solve. (Re-)set velocity boundary conditions just to be on the safe s...
double Re_invFr
Reynolds/Froude number.
void doc_solution(DocInfo &doc_info)
Doc the solution.
void actions_after_adapt()
After adaptation: Unpin pressure and pin redudant pressure dofs.
void zero_body_force(const double &time, const Vector< double > &x, Vector< double > &result)
Zero functional body force.
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.
~QuarterCircleDrivenCavityProblem2()
Destructor: Empty.