full_tube.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: 1207 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-06-24 15:10:32 +0100 (Fri, 24 Jun 2016) $
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 steady entry flow problem in a
31 ///complete straight tube
32 
33 //Generic routines
34 #include "generic.h"
35 #include "navier_stokes.h"
36 
37 // The mesh
38 #include "meshes/tube_mesh.h"
39 
40 using namespace std;
41 
42 using namespace oomph;
43 
44 
45 //=start_of_MyCylinder===================================
46 //A geometric object that represents the geometry of the domain
47 //================================================================
48 class MyCylinder : public GeomObject
49 {
50 public:
51 
52  /// Constructor that takes the radius of the tube as its argument
53  MyCylinder(const double &radius) :
54  GeomObject(3,3), Radius(radius) { }
55 
56 /// Destructor
57 virtual~MyCylinder(){}
58 
59 ///Lagrangian coordinate xi
60 void position (const Vector<double>& xi, Vector<double>& r) const
61 {
62  r[0] = xi[2]*Radius*cos(xi[1]);
63  r[1] = xi[0];
64  r[2] = -xi[2]*Radius*sin(xi[1]);
65 }
66 
67 
68 /// Return the position of the tube as a function of time
69 /// (doesn't move as a function of time)
70 void position(const unsigned& t,
71  const Vector<double>& xi, Vector<double>& r) const
72  {
73  position(xi,r);
74  }
75 
76 private:
77 
78  ///Storage for the radius of the tube
79  double Radius;
80 
81 };
82 
83 //=start_of_namespace================================================
84 /// Namespace for physical parameters
85 //===================================================================
86 namespace Global_Physical_Variables
87 {
88  /// Reynolds number
89  double Re=25;
90 
91 } // end_of_namespace
92 
93 
94 //=start_of_problem_class=============================================
95 /// Entry flow problem in tapered tube domain
96 //====================================================================
97 template<class ELEMENT>
98 class SteadyTubeProblem : public Problem
99 {
100 
101 public:
102 
103  /// Constructor: Pass DocInfo object and target errors
104  SteadyTubeProblem(DocInfo& doc_info, const double& min_error_target,
105  const double& max_error_target);
106 
107  /// Destructor (empty)
109 
110  /// \short Update the problem specs before solve
111  void actions_before_newton_solve();
112 
113  /// After adaptation: Pin redudant pressure dofs.
115  {
116  // Pin redudant pressure dofs
117  RefineableNavierStokesEquations<3>::
118  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
119  }
120 
121  /// Doc the solution
122  void doc_solution();
123 
124  /// \short Overload generic access function by one that returns
125  /// a pointer to the specific mesh
126  RefineableTubeMesh<ELEMENT>* mesh_pt()
127  {
128  return dynamic_cast<RefineableTubeMesh<ELEMENT>*>(Problem::mesh_pt());
129  }
130 
131 private:
132 
133  /// Doc info object
134  DocInfo Doc_info;
135 
136  ///Pointer to GeomObject that specifies the domain volume
137  GeomObject *Volume_pt;
138 
139 }; // end_of_problem_class
140 
141 
142 
143 
144 //=start_of_constructor===================================================
145 /// Constructor: Pass DocInfo object and error targets
146 //========================================================================
147 template<class ELEMENT>
149  const double& min_error_target,
150  const double& max_error_target)
151  : Doc_info(doc_info)
152 {
153 
154  // Setup mesh:
155  //------------
156 
157  // Create GeomObject that specifies the domain geometry
158  //The radius of the tube is one
159  Volume_pt=new MyCylinder(1.0);
160 
161  //Define pi
162  const double pi = MathematicalConstants::Pi;
163 
164  //Set the centerline coordinates spanning the mesh
165  //Tube is on length pi
166  Vector<double> centreline_limits(2);
167  centreline_limits[0] = 0.0;
168  centreline_limits[1] = pi;
169 
170  //Set the positions of the angles that divide the outer ring
171  //These must be in the range -pi,pi, ordered from smallest to
172  //largest
173  Vector<double> theta_positions(4);
174  theta_positions[0] = -0.75*pi;
175  theta_positions[1] = -0.25*pi;
176  theta_positions[2] = 0.25*pi;
177  theta_positions[3] = 0.75*pi;
178 
179  //Define the radial fraction of the central box (always halfway
180  //along the radius)
181  Vector<double> radial_frac(4,0.5);
182 
183  // Number of layers in the initial mesh
184  unsigned nlayer=6;
185 
186  // Build and assign mesh
187  Problem::mesh_pt()= new RefineableTubeMesh<ELEMENT>(Volume_pt,
188  centreline_limits,
189  theta_positions,
190  radial_frac,
191  nlayer);
192 
193  // Set error estimator
194  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
195  mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
196 
197  // Error targets for adaptive refinement
198  mesh_pt()->max_permitted_error()=max_error_target;
199  mesh_pt()->min_permitted_error()=min_error_target;
200 
201  // Set the boundary conditions for this problem: All nodal values are
202  // free by default -- just pin the ones that have Dirichlet conditions
203  // here.
204  //Choose the conventional form by setting gamma to zero
205  //The boundary conditions will be pseudo-traction free (d/dn = 0)
206  ELEMENT::Gamma[0] = 0.0;
207  ELEMENT::Gamma[1] = 0.0;
208  ELEMENT::Gamma[2] = 0.0;
209 
210  //Loop over the boundaries
211  unsigned num_bound = mesh_pt()->nboundary();
212  for(unsigned ibound=0;ibound<num_bound;ibound++)
213  {
214  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
215  for (unsigned inod=0;inod<num_nod;inod++)
216  {
217  // Boundary 0 is the inlet symmetry boundary:
218  // Boundary 1 is the tube wall
219  // Pin all values
220  if((ibound==0) || (ibound==1))
221  {
222  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
223  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
224  mesh_pt()->boundary_node_pt(ibound,inod)->pin(2);
225  }
226  }
227  } // end loop over boundaries
228 
229  // Loop over the elements to set up element-specific
230  // things that cannot be handled by constructor
231  unsigned n_element = mesh_pt()->nelement();
232  for(unsigned i=0;i<n_element;i++)
233  {
234  // Upcast from GeneralisedElement to the present element
235  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
236 
237  //Set the Reynolds number, etc
238  el_pt->re_pt() = &Global_Physical_Variables::Re;
239  }
240 
241  // Pin redudant pressure dofs
242  RefineableNavierStokesEquations<3>::
243  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
244 
245  //Attach the boundary conditions to the mesh
246  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
247 
248 } // end_of_constructor
249 
250 
251 //=start_of_actions_before_newton_solve===================================
252 /// Set the inflow boundary conditions
253 //========================================================================
254 template<class ELEMENT>
256 {
257 
258  // (Re-)assign velocity profile at inflow values
259  //--------------------------------------------
260  unsigned ibound=0;
261  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
262  for (unsigned inod=0;inod<num_nod;inod++)
263  {
264  // Recover coordinates of tube relative to centre position
265  double x=mesh_pt()->boundary_node_pt(ibound,inod)->x(0);
266  double z=mesh_pt()->boundary_node_pt(ibound,inod)->x(2);
267  //Calculate the radius
268  double r=sqrt(x*x+z*z);
269 
270  // Blunt profile for axial velocity (component 1 at the inlet)
271  mesh_pt()->boundary_node_pt(ibound,inod)->
272  set_value(1,(1.0-pow(r,10.0)));
273  }
274 
275 } // end_of_actions_before_newton_solve
276 
277 
278 //=start_of_doc_solution==================================================
279 /// Doc the solution
280 //========================================================================
281 template<class ELEMENT>
283 {
284  //Output file stream
285  ofstream some_file;
286  char filename[100];
287 
288  // Number of plot points
289  unsigned npts;
290  npts=5;
291 
292  //Need high precision for large radii of curvature
293  //some_file.precision(10);
294  // Output solution labelled by the Reynolds number
295  sprintf(filename,"%s/soln_Re%g.dat",Doc_info.directory().c_str(),
297  some_file.open(filename);
298  mesh_pt()->output(some_file,npts);
299  some_file.close();
300 
301 } // end_of_doc_solution
302 
303 
304 
305 
306 ////////////////////////////////////////////////////////////////////////
307 ////////////////////////////////////////////////////////////////////////
308 ////////////////////////////////////////////////////////////////////////
309 
310 
311 //=start_of_main=======================================================
312 /// Driver for 3D entry flow into a straight tube.
313 //=====================================================================
314 int main(int argc, char* argv[])
315 {
316  //The purpose of this code is to test that the position of nodes on
317  //curvilinear boundaries are correctly updated after nodes that were
318  //hanging become non-hanging.
319  unsigned max_adapt=1;
320  double max_error_target=0.02 , min_error_target=0.002;
321 
322 
323  // Set up doc info
324  DocInfo doc_info;
325 
326  // Do Taylor-Hood elements
327  //------------------------
328  {
329  // Set output directory
330  doc_info.set_directory("RESLT_TH");
331 
332  // Step number
333  doc_info.number()=0;
334 
335  // Build problem
337  problem(doc_info,min_error_target,max_error_target);
338 
339  cout << " Doing Taylor-Hood elements " << std::endl;
340 
341  for(unsigned i=0;i<2;i++)
342  {
343  // Solve the problem
344  problem.newton_solve(max_adapt);
345  // Doc solution after solving
346  problem.doc_solution();
347 
349  ++doc_info.number();
350  }
351  }
352 
353  exit(1);
354  // Do Crouzeix-Raviart elements
355  //------------------------
356  {
357  // Set output directory
358  doc_info.set_directory("RESLT_CR");
359 
360  // Step number
361  doc_info.number()=0;
362 
363  // Build problem
365  problem(doc_info,min_error_target,max_error_target);
366 
367  cout << " Doing Crouzeix-Raviart elements " << std::endl;
368 
369  for(unsigned i=0;i<2;i++)
370  {
371  // Solve the problem
372  problem.newton_solve(max_adapt);
373  // Doc solution after solving
374  problem.doc_solution();
375 
377  ++doc_info.number();
378  }
379  }
380 } // end_of_main
381 
382 
int main(int argc, char *argv[])
Driver for 3D entry flow into a straight tube.
Definition: full_tube.cc:314
~SteadyTubeProblem()
Destructor (empty)
Definition: full_tube.cc:108
MyCylinder(const double &radius)
Constructor that takes the radius of the tube as its argument.
Definition: full_tube.cc:53
double Radius
Storage for the radius of the tube.
Definition: full_tube.cc:79
virtual ~MyCylinder()
Destructor.
Definition: full_tube.cc:57
Entry flow problem in tapered tube domain.
Definition: full_tube.cc:98
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Definition: full_tube.cc:70
SteadyTubeProblem(DocInfo &doc_info, const double &min_error_target, const double &max_error_target)
Constructor: Pass DocInfo object and target errors.
Definition: full_tube.cc:148
double Re
Reynolds number.
Definition: full_tube.cc:89
void doc_solution()
Doc the solution.
Definition: full_tube.cc:282
GeomObject * Volume_pt
Pointer to GeomObject that specifies the domain volume.
Definition: full_tube.cc:137
RefineableTubeMesh< ELEMENT > * mesh_pt()
Overload generic access function by one that returns a pointer to the specific mesh.
Definition: full_tube.cc:126
void actions_after_adapt()
After adaptation: Pin redudant pressure dofs.
Definition: full_tube.cc:114
void actions_before_newton_solve()
Update the problem specs before solve.
Definition: full_tube.cc:255
void position(const Vector< double > &xi, Vector< double > &r) const
Lagrangian coordinate xi.
Definition: full_tube.cc:60
DocInfo Doc_info
Doc info object.
Definition: full_tube.cc:134