three_d_entry_flow.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 a 3D navier stokes entry flow problem in quarter tube domain
31 
32 //Generic routines
33 #include "generic.h"
34 #include "navier_stokes.h"
35 
36 // The mesh
37 #include "meshes/quarter_tube_mesh.h"
38 
39 using namespace std;
40 
41 using namespace oomph;
42 
43 //=start_of_namespace================================================
44 /// Namespace for physical parameters
45 //===================================================================
46 namespace Global_Physical_Variables
47 {
48  /// Reynolds number
49  double Re=100;
50 } // end_of_namespace
51 
52 
53 //=start_of_problem_class=============================================
54 /// Entry flow problem in quarter tube domain
55 //====================================================================
56 template<class ELEMENT>
57 class EntryFlowProblem : public Problem
58 {
59 
60 public:
61 
62  /// Constructor: Pass DocInfo object and target errors
63  EntryFlowProblem(DocInfo& doc_info, const double& min_error_target,
64  const double& max_error_target);
65 
66  /// Destructor (empty)
68 
69  /// Doc the solution after solve
71  {
72  // Doc solution after solving
73  doc_solution();
74 
75  // Increment label for output files
76  Doc_info.number()++;
77  }
78 
79  /// \short Update the problem specs before solve
80  void actions_before_newton_solve();
81 
82  /// After adaptation: Pin redudant pressure dofs.
84  {
85  // Pin redudant pressure dofs
86  RefineableNavierStokesEquations<3>::
87  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
88  }
89 
90  /// Doc the solution
91  void doc_solution();
92 
93  /// \short Overload generic access function by one that returns
94  /// a pointer to the specific mesh
95  RefineableQuarterTubeMesh<ELEMENT>* mesh_pt()
96  {
97  return dynamic_cast<RefineableQuarterTubeMesh<ELEMENT>*>(Problem::mesh_pt());
98  }
99 
100 private:
101 
102  /// Exponent for bluntness of velocity profile
103  int Alpha;
104 
105  /// Doc info object
106  DocInfo Doc_info;
107 
108 }; // end_of_problem_class
109 
110 
111 
112 
113 //=start_of_constructor===================================================
114 /// Constructor: Pass DocInfo object and error targets
115 //========================================================================
116 template<class ELEMENT>
118  const double& min_error_target,
119  const double& max_error_target)
120  : Doc_info(doc_info)
121 {
122 
123  // Setup mesh:
124  //------------
125 
126  // Create geometric objects: Elliptical tube with half axes = radius = 1.0
127  double radius=1.0;
128  GeomObject* Wall_pt=new EllipticalTube(radius,radius);
129 
130  // Boundaries on object
131  Vector<double> xi_lo(2);
132  // height of inflow
133  xi_lo[0]=0.0;
134  // start of Wall_pt
135  xi_lo[1]=0.0;
136 
137  Vector<double> xi_hi(2);
138  // height of outflow
139  xi_hi[0]=7.0;
140  // end of Wall_pt
141  xi_hi[1]=2.0*atan(1.0);
142 
143  // # of layers
144  unsigned nlayer=6;
145 
146  //Radial divider is located half-way along the circumference
147  double frac_mid=0.5;
148 
149  // Build and assign mesh
150  Problem::mesh_pt()=
151  new RefineableQuarterTubeMesh<ELEMENT>(Wall_pt,xi_lo,frac_mid,xi_hi,nlayer);
152 
153 
154  // Set error estimator
155  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
156  mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
157 
158  // Error targets for adaptive refinement
159  mesh_pt()->max_permitted_error()=max_error_target;
160  mesh_pt()->min_permitted_error()=min_error_target;
161 
162 
163 
164  //Doc the boundaries
165  ofstream some_file;
166  char filename[100];
167  sprintf(filename,"boundaries.dat");
168  some_file.open(filename);
169  mesh_pt()->output_boundaries(some_file);
170  some_file.close();
171 
172 
173  // Set the boundary conditions for this problem: All nodal values are
174  // free by default -- just pin the ones that have Dirichlet conditions
175  // here.
176  unsigned num_bound = mesh_pt()->nboundary();
177  for(unsigned ibound=0;ibound<num_bound;ibound++)
178  {
179  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
180  for (unsigned inod=0;inod<num_nod;inod++)
181  {
182  // Boundary 1 is the vertical symmetry boundary: We allow flow in
183  //the y-direction. Elsewhere, pin the y velocity
184  if(ibound!=1) mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
185 
186  // Boundary 2 is the horizontal symmetry boundary: We allow flow in
187  //the x-direction. Elsewhere, pin the x velocity
188  if(ibound!=2) mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
189 
190  // Boundaries 0 and 3 are the inflow and the wall respectively.
191  // Pin the axial velocity because of the prescribed inflow
192  // profile and the no slip on the stationary wall, respectively
193  if((ibound==0) || (ibound==3))
194  {
195  mesh_pt()->boundary_node_pt(ibound,inod)->pin(2);
196  }
197  }
198  } // end loop over boundaries
199 
200  // Loop over the elements to set up element-specific
201  // things that cannot be handled by constructor
202  unsigned n_element = mesh_pt()->nelement();
203  for(unsigned i=0;i<n_element;i++)
204  {
205  // Upcast from GeneralisedElement to the present element
206  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
207 
208  //Set the Reynolds number, etc
209  el_pt->re_pt() = &Global_Physical_Variables::Re;
210  }
211 
212  // Pin redudant pressure dofs
213  RefineableNavierStokesEquations<3>::
214  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
215 
216  // Set Poiseuille flow as initial for solution
217  // Inflow will be overwritten in actions_before_newton_solve()
218  unsigned n_nod=mesh_pt()->nnode();
219  for (unsigned j=0;j<n_nod;j++)
220  {
221  Node* node_pt=mesh_pt()->node_pt(j);
222  // Recover coordinates
223  double x=node_pt->x(0);
224  double y=node_pt->x(1);
225  double r=sqrt(x*x+y*y );
226 
227  // Poiseuille flow
228  node_pt->set_value(0,0.0);
229  node_pt->set_value(1,0.0);
230  node_pt->set_value(2,(1.0-r*r));
231  }
232 
233  // Set the exponent for bluntness: Alpha=2 --> Poisseuille; anything
234  // larger makes the inflow blunter
235  Alpha=20;
236 
237  //Attach the boundary conditions to the mesh
238  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
239 
240 } // end_of_constructor
241 
242 
243 //=start_of_actions_before_newton_solve==========================================
244 /// Set the inflow boundary conditions
245 //========================================================================
246 template<class ELEMENT>
248 {
249 
250  // (Re-)assign velocity profile at inflow values
251  //--------------------------------------------
252 
253  // Setup bluntish parallel inflow on boundary 0:
254  unsigned ibound=0;
255  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
256  for (unsigned inod=0;inod<num_nod;inod++)
257  {
258  // Recover coordinates
259  double x=mesh_pt()->boundary_node_pt(ibound,inod)->x(0);
260  double y=mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
261  double r=sqrt(x*x+y*y);
262 
263  // Bluntish profile for axial velocity (component 2)
264  mesh_pt()->boundary_node_pt(ibound,inod)->
265  set_value(2,(1.0-pow(r,Alpha)));
266  }
267 
268 } // end_of_actions_before_newton_solve
269 
270 
271 //=start_of_doc_solution==================================================
272 /// Doc the solution
273 //========================================================================
274 template<class ELEMENT>
276 {
277 
278  ofstream some_file;
279  char filename[100];
280 
281  // Number of plot points
282  unsigned npts;
283  npts=5;
284 
285  // Output solution
286  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
287  Doc_info.number());
288  some_file.open(filename);
289  mesh_pt()->output(some_file,npts);
290  some_file.close();
291 
292 } // end_of_doc_solution
293 
294 
295 
296 
297 ////////////////////////////////////////////////////////////////////////
298 ////////////////////////////////////////////////////////////////////////
299 ////////////////////////////////////////////////////////////////////////
300 
301 
302 //=start_of_main=======================================================
303 /// Driver for 3D entry flow into a quarter tube. If there are
304 /// any command line arguments, we regard this as a validation run
305 /// and perform only a single adaptation
306 //=====================================================================
307 int main(int argc, char* argv[])
308 {
309 
310  // Store command line arguments
311  CommandLineArgs::setup(argc,argv);
312 
313  // Allow (up to) five rounds of fully automatic adapation in response to
314  //-----------------------------------------------------------------------
315  // error estimate
316  //---------------
317  unsigned max_adapt;
318  double max_error_target,min_error_target;
319 
320  // Set max number of adaptations in black-box Newton solver and
321  // error targets for adaptation
322  if (CommandLineArgs::Argc==1)
323  {
324  // Up to five adaptations
325  max_adapt=5;
326 
327  // Error targets for adaptive refinement
328  max_error_target=0.005;
329  min_error_target=0.0005;
330  }
331  // Validation run: Only one adaptation. Relax error targets
332  // to ensure that not all elements are refined so we're getting
333  // some hanging nodes.
334  else
335  {
336  // Validation run: Just one round of adaptation
337  max_adapt=1;
338 
339  // Error targets for adaptive refinement
340  max_error_target=0.02;
341  min_error_target=0.002;
342  }
343  // end max_adapt setup
344 
345 
346  // Set up doc info
347  DocInfo doc_info;
348 
349  // Do Taylor-Hood elements
350  //------------------------
351  {
352  // Set output directory
353  doc_info.set_directory("RESLT_TH");
354 
355  // Step number
356  doc_info.number()=0;
357 
358  // Build problem
360  problem(doc_info,min_error_target,max_error_target);
361 
362  cout << " Doing Taylor-Hood elements " << std::endl;
363 
364 
365  // Doc solution after solving
366  problem.doc_solution();
367 
368  // Solve the problem
369  problem.newton_solve(max_adapt);
370  }
371 
372 
373  // Do Crouzeix-Raviart elements
374  //-----------------------------
375  {
376  // Set output directory
377  doc_info.set_directory("RESLT_CR");
378 
379  // Step number
380  doc_info.number()=0;
381 
382  // Build problem
384  problem(doc_info,min_error_target,max_error_target);
385 
386  cout << " Doing Crouzeix-Raviart elements " << std::endl;
387 
388  // Solve the problem
389  problem.newton_solve(max_adapt);
390  }
391 
392 } // end_of_main
393 
394 
int main(int argc, char *argv[])
void actions_after_newton_solve()
Doc the solution after solve.
EntryFlowProblem(DocInfo &doc_info, const double &min_error_target, const double &max_error_target)
Constructor: Pass DocInfo object and target errors.
void actions_after_adapt()
After adaptation: Pin redudant pressure dofs.
Entry flow problem in quarter tube domain.
int Alpha
Exponent for bluntness of velocity profile.
~EntryFlowProblem()
Destructor (empty)
double Re
Reynolds number.
Definition: full_tube.cc:89
DocInfo Doc_info
Doc info object.
void actions_before_newton_solve()
Update the problem specs before solve.
RefineableQuarterTubeMesh< ELEMENT > * mesh_pt()
Overload generic access function by one that returns a pointer to the specific mesh.
void doc_solution()
Doc the solution.