rayleigh_traction_channel.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: 1292 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-03-06 18:57:19 +0000 (Mon, 06 Mar 2017) $
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 "Rayleigh channel" driven by an applied traction
31 
32 // The oomphlib headers
33 #include "generic.h"
34 #include "navier_stokes.h"
35 
36 // The mesh
37 #include "meshes/rectangular_quadmesh.h"
38 
39 using namespace std;
40 
41 using namespace oomph;
42 
43 //===start_of_namespace=================================================
44 /// Namepspace for global parameters
45 //======================================================================
46 namespace Global_Parameters
47 {
48  /// Reynolds number
49  double Re;
50 
51  /// Womersley = Reynolds times Strouhal
52  double ReSt;
53 
54  /// Flag for long/short run: Default = perform long run
55  unsigned Long_run_flag=1;
56 
57  /// \short Flag for impulsive start: Default = start from exact
58  /// time-periodic solution.
60 
61 } // end of namespace
62 
63 
64 //==start_of_exact_solution=============================================
65 /// Namespace for exact solution
66 //======================================================================
67 namespace ExactSoln
68 {
69 
70  /// Exact solution of the problem as a vector
71  void get_exact_u(const double& t, const Vector<double>& x, Vector<double>& u)
72  {
73  double y=x[1];
74  // I=sqrt(-1)
75  complex<double> I(0.0,1.0);
76  // omega
77  double omega=2.0*MathematicalConstants::Pi;
78  // lambda
79  complex<double> lambda(0.0,omega*Global_Parameters::ReSt);
80  lambda = I*sqrt(lambda);
81 
82  // Solution
83  complex<double> sol(
84  exp(complex<double>(0.0,omega*t)) *
85  (exp(lambda*complex<double>(0.0,y))-exp(lambda*complex<double>(0.0,-y)))
86  /(exp(I*lambda)-exp(-I*lambda)) );
87 
88  // Extract real part and stick into vector
89  u.resize(2);
90  u[0]=real(sol);
91  u[1]=0.0;
92 
93  }
94 
95  /// Traction required at the top boundary
96  void prescribed_traction(const double& t,
97  const Vector<double>& x,
98  const Vector<double> &n,
99  Vector<double>& traction)
100  {
101  double y=x[1];
102  // I=sqrt(-1)
103  complex<double> I(0.0,1.0);
104  // omega
105  double omega=2.0*MathematicalConstants::Pi;
106  // lambda
107  complex<double> lambda(0.0,omega*Global_Parameters::ReSt);
108  lambda = I*sqrt(lambda);
109 
110  // Solution
111  complex<double> sol(
112  exp(complex<double>(0.0,omega*t)) *
113  (exp(lambda*complex<double>(0.0,y))+exp(lambda*complex<double>(0.0,-y)))
114  *I*lambda/(exp(I*lambda)-exp(-I*lambda)) );
115 
116  //Extract real part and stick into vector
117  traction.resize(2);
118  traction[0]=real(sol);
119  traction[1]=0.0;
120  }
121 
122 } // end of exact_solution
123 
124 
125 //===start_of_problem_class=============================================
126 /// Rayleigh-type problem: 2D channel flow driven by traction bc
127 //======================================================================
128 template<class ELEMENT, class TIMESTEPPER>
129 class RayleighTractionProblem : public Problem
130 {
131 public:
132 
133  /// Constructor: Pass number of elements in x and y directions and
134  /// lengths
135  RayleighTractionProblem(const unsigned &nx, const unsigned &ny,
136  const double &lx, const double &ly);
137 
138  /// Update before solve is empty
140 
141  /// \short Update after solve is empty
143 
144  /// Actions before timestep (empty)
146 
147  /// Run an unsteady simulation
148  void unsteady_run(DocInfo& doc_info);
149 
150  /// Doc the solution
151  void doc_solution(DocInfo& doc_info);
152 
153  /// \short Set initial condition (incl previous timesteps) according
154  /// to specified function.
155  void set_initial_condition();
156 
157  /// \short Create traction elements on boundary b of the Mesh pointed
158  /// to by bulk_mesh_pt and add them to the Mesh object pointed to by
159  /// surface_mesh_pt
160  void create_traction_elements(const unsigned &b,
161  Mesh* const &bulk_mesh_pt,
162  Mesh* const &surface_mesh_pt);
163 
164 private:
165 
166  /// Pointer to the "bulk" mesh
167  RectangularQuadMesh<ELEMENT>* Bulk_mesh_pt;
168 
169  /// Pointer to the "surface" mesh
171 
172  /// Trace file
173  ofstream Trace_file;
174 
175 }; // end of problem_class
176 
177 
178 
179 //===start_of_constructor=============================================
180 /// Problem constructor
181 //====================================================================
182 template<class ELEMENT,class TIMESTEPPER>
184 (const unsigned &nx, const unsigned &ny,
185  const double &lx, const double& ly)
186 {
187  //Allocate the timestepper
188  add_time_stepper_pt(new TIMESTEPPER);
189 
190  //Now create the mesh with periodic boundary conditions in x direction
191  bool periodic_in_x=true;
192  Bulk_mesh_pt =
193  new RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,periodic_in_x,
194  time_stepper_pt());
195 
196  // Create "surface mesh" that will contain only the prescribed-traction
197  // elements. The constructor just creates the mesh without
198  // giving it any elements, nodes, etc.
199  Surface_mesh_pt = new Mesh;
200 
201  // Create prescribed-traction elements from all elements that are
202  // adjacent to boundary 2, but add them to a separate mesh.
203  create_traction_elements(2,Bulk_mesh_pt,Surface_mesh_pt);
204 
205  // Add the two sub meshes to the problem
206  add_sub_mesh(Bulk_mesh_pt);
207  add_sub_mesh(Surface_mesh_pt);
208 
209  // Combine all submeshes into a single Mesh
210  build_global_mesh();
211 
212  // Set the boundary conditions for this problem: All nodes are
213  // free by default -- just pin the ones that have Dirichlet conditions
214  // here
215  unsigned num_bound=Bulk_mesh_pt->nboundary();
216  for(unsigned ibound=0;ibound<num_bound;ibound++)
217  {
218  unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
219  for (unsigned inod=0;inod<num_nod;inod++)
220  {
221  // No slip on bottom
222  // (DO NOT PIN TOP BOUNDARY, SINCE TRACTION DEFINES ITS VELOCITY!)
223  if (ibound==0)
224  {
225  Bulk_mesh_pt->boundary_node_pt(ibound,inod)->pin(0);
226  Bulk_mesh_pt->boundary_node_pt(ibound,inod)->pin(1);
227  }
228  // Horizontal outflow on the left (and right -- right bc not
229  // strictly necessary because of symmetry)
230  else if ((ibound==1) || (ibound==3))
231  {
232  Bulk_mesh_pt->boundary_node_pt(ibound,inod)->pin(1);
233  }
234  }
235  } // end loop over boundaries
236 
237  //Complete the problem setup to make the elements fully functional
238 
239  //Loop over the elements
240  unsigned n_el = Bulk_mesh_pt->nelement();
241  for(unsigned e=0;e<n_el;e++)
242  {
243  //Cast to a fluid element
244  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
245 
246  //Set the Reynolds number, etc
247  el_pt->re_pt() = &Global_Parameters::Re;
248  el_pt->re_st_pt() = &Global_Parameters::ReSt;
249  }
250 
251  // Loop over the flux elements to pass pointer to prescribed traction function
252  // and pointer to global time object
253  n_el=Surface_mesh_pt->nelement();
254  for(unsigned e=0;e<n_el;e++)
255  {
256  // Upcast from GeneralisedElement to traction element
257  NavierStokesTractionElement<ELEMENT> *el_pt =
258  dynamic_cast< NavierStokesTractionElement<ELEMENT>*>(
259  Surface_mesh_pt->element_pt(e));
260 
261  // Set the pointer to the prescribed traction function
262  el_pt->traction_fct_pt() =
264  }
265 
266  //Assgn equation numbers
267  cout << assign_eqn_numbers() << std::endl;
268 
269 } // end of constructor
270 
271 //======================start_of_set_initial_condition====================
272 /// \short Set initial condition: Assign previous and current values
273 /// from exact solution.
274 //========================================================================
275 template<class ELEMENT,class TIMESTEPPER>
277 {
278  // Backup time in global Time object
279  double backed_up_time=time_pt()->time();
280 
281  // Past history needs to be established for t=time0-deltat, ...
282  // Then provide current values (at t=time0) which will also form
283  // the initial guess for the first solve at t=time0+deltat
284 
285  // Vector of exact solution value
286  Vector<double> soln(2);
287  Vector<double> x(2);
288 
289  //Find number of nodes in mesh
290  unsigned num_nod = mesh_pt()->nnode();
291 
292  // Set continuous times at previous timesteps:
293  // How many previous timesteps does the timestepper use?
294  int nprev_steps=time_stepper_pt()->nprev_values();
295  Vector<double> prev_time(nprev_steps+1);
296  for (int t=nprev_steps;t>=0;t--)
297  {
298  prev_time[t]=time_pt()->time(unsigned(t));
299  }
300 
301  // Loop over current & previous timesteps
302  for (int t=nprev_steps;t>=0;t--)
303  {
304  // Continuous time
305  double time=prev_time[t];
306  cout << "setting IC at time =" << time << std::endl;
307 
308  // Loop over the nodes to set initial guess everywhere
309  for (unsigned n=0;n<num_nod;n++)
310  {
311  // Get nodal coordinates
312  x[0]=mesh_pt()->node_pt(n)->x(0);
313  x[1]=mesh_pt()->node_pt(n)->x(1);
314 
315  // Get exact solution at previous time
316  ExactSoln::get_exact_u(time,x,soln);
317 
318  // Assign solution
319  mesh_pt()->node_pt(n)->set_value(t,0,soln[0]);
320  mesh_pt()->node_pt(n)->set_value(t,1,soln[1]);
321 
322  // Loop over coordinate directions: Mesh doesn't move, so
323  // previous position = present position
324  for (unsigned i=0;i<2;i++)
325  {
326  mesh_pt()->node_pt(n)->x(t,i)=x[i];
327  }
328  }
329  }
330 
331  // Reset backed up time for global timestepper
332  time_pt()->time()=backed_up_time;
333 
334 } // end of set_initial_condition
335 
336 
337 
338 
339 //==start_of_doc_solution=================================================
340 /// Doc the solution
341 //========================================================================
342 template<class ELEMENT,class TIMESTEPPER>
344 {
345 
346  ofstream some_file;
347  char filename[100];
348 
349  // Number of plot points
350  unsigned npts=5;
351 
352 
353  // Output solution
354  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
355  doc_info.number());
356  some_file.open(filename);
357  Bulk_mesh_pt->output(some_file,npts);
358 
359  // Write file as a tecplot text object
360  some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
361  << time_pt()->time() << "\"";
362  // ...and draw a horizontal line whose length is proportional
363  // to the elapsed time
364  some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
365  some_file << "1" << std::endl;
366  some_file << "2" << std::endl;
367  some_file << " 0 0" << std::endl;
368  some_file << time_pt()->time()*20.0 << " 0" << std::endl;
369 
370  some_file.close();
371 
372  // Output exact solution
373  //----------------------
374  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
375  doc_info.number());
376  some_file.open(filename);
377  Bulk_mesh_pt->output_fct(some_file,npts,time_pt()->time(),
379  some_file.close();
380 
381  // Doc error
382  //----------
383  double error,norm;
384  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
385  doc_info.number());
386  some_file.open(filename);
387  Bulk_mesh_pt->compute_error(some_file,
389  time_pt()->time(),
390  error,norm);
391  some_file.close();
392 
393  // Doc solution and error
394  //-----------------------
395  cout << "error: " << error << std::endl;
396  cout << "norm : " << norm << std::endl << std::endl;
397 
398  // Get time, position and exact soln
399  Vector<double> x(2), u(2);
400  double time=time_pt()->time();
401  Node* node_pt=dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(37))->node_pt(1);
402  x[0] = node_pt->x(0);
403  x[1] = node_pt->x(1);
404  ExactSoln::get_exact_u(time,x,u);
405 
406  // doc
407  Trace_file << time << " "
408  << x[0] << " "
409  << x[1] << " "
410  << node_pt->value(0) << " "
411  << node_pt->value(1) << " "
412  << u[0] << " "
413  << u[1] << " "
414  << abs(u[0]-node_pt->value(0)) << " "
415  << abs(u[1]-node_pt->value(1)) << " "
416  << error << " "
417  << norm << " "
418  << std::endl;
419 
420 } // end_of_doc_solution
421 
422 //============start_of_create_traction_elements==========================
423 /// Create Navier-Stokes traction elements on the b-th boundary of the
424 /// Mesh object pointed to by bulk_mesh_pt and add the elements to the
425 /// Mesh object pointeed to by surface_mesh_pt.
426 //=======================================================================
427 template<class ELEMENT,class TIMESTEPPER>
429 create_traction_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
430  Mesh* const &surface_mesh_pt)
431 {
432  // How many bulk elements are adjacent to boundary b?
433  unsigned n_element = bulk_mesh_pt->nboundary_element(b);
434 
435  // Loop over the bulk elements adjacent to boundary b?
436  for(unsigned e=0;e<n_element;e++)
437  {
438  // Get pointer to the bulk element that is adjacent to boundary b
439  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
440  bulk_mesh_pt->boundary_element_pt(b,e));
441 
442  //What is the index of the face of element e along boundary b
443  int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
444 
445  // Build the corresponding prescribed-flux element
446  NavierStokesTractionElement<ELEMENT>* flux_element_pt = new
447  NavierStokesTractionElement<ELEMENT>(bulk_elem_pt,face_index);
448 
449  //Add the prescribed-flux element to the surface mesh
450  surface_mesh_pt->add_element_pt(flux_element_pt);
451 
452  } //end of loop over bulk elements adjacent to boundary b
453 
454 } // end of create_traction_elements
455 
456 
457 
458 //===start_of_unsteady_run=====================================================
459 /// Unsteady run...
460 //=============================================================================
461 template<class ELEMENT,class TIMESTEPPER>
463 {
464 
465  // Open trace file
466  char filename[100];
467  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
468  Trace_file.open(filename);
469 
470  // initialise Trace file
471  Trace_file << "time" << ", "
472  << "x" << ", "
473  << "y" << ", "
474  << "u_1" << ", "
475  << "u_2" << ", "
476  << "u_exact_1" << ", "
477  << "u_exact_2" << ", "
478  << "error_1" << ", "
479  << "error_2" << ", "
480  << "L2 error" << ", "
481  << "L2 norm" << ", " << std::endl;
482 
483  //Set value of dt
484  double dt = 0.025;
485 
487  {
488  // Initialise all history values for an impulsive start
489  assign_initial_values_impulsive(dt);
490  cout << "IC = impulsive start" << std::endl;
491  }
492  else
493  {
494  // Initialise timestep
495  initialise_dt(dt);
496  // Set initial conditions.
497  set_initial_condition();
498  cout << "IC = exact solution" << std::endl;
499  }
500 
501  //Now do many timesteps
502  unsigned ntsteps=500;
503 
504  // If validation run only do 5 timesteps
506  {
507  ntsteps=5;
508  cout << "validation run" << std::endl;
509  }
510 
511  // Doc initial condition
512  doc_solution(doc_info);
513 
514  // increment counter
515  doc_info.number()++;
516 
517  //Loop over the timesteps
518  for(unsigned t=1;t<=ntsteps;t++)
519  {
520  cout << "TIMESTEP " << t << std::endl;
521 
522  //Take one fixed timestep
523  unsteady_newton_solve(dt);
524 
525  //Output the time
526  cout << "Time is now " << time_pt()->time() << std::endl;
527 
528  // Doc solution
529  doc_solution(doc_info);
530 
531  // increment counter
532  doc_info.number()++;
533  }
534 
535 } // end of unsteady run
536 
537 
538 //===start_of_main======================================================
539 /// Driver code for Rayleigh-type problem
540 //======================================================================
541 int main(int argc, char* argv[])
542 {
543 
544 
545  /// Convert command line arguments (if any) into flags:
546  if (argc==1)
547  {
548  cout << "No command line arguments specified -- using defaults." << std::endl;
549  }
550  else if (argc==3)
551  {
552  cout << "Two command line arguments specified:" << std::endl;
553  // Flag for long run
554  Global_Parameters::Long_run_flag=atoi(argv[1]);
555  /// Flag for impulsive start
557  }
558  else
559  {
560  std::string error_message =
561  "Wrong number of command line arguments. Specify none or two.\n";
562  error_message +=
563  "Arg1: Long_run_flag [0/1]\n";
564  error_message +=
565  "Arg2: Impulsive_start_flag [0/1]\n";
566 
567  throw OomphLibError(error_message,
568  OOMPH_CURRENT_FUNCTION,
569  OOMPH_EXCEPTION_LOCATION);
570  }
571  cout << "Long run flag: "
572  << Global_Parameters::Long_run_flag << std::endl;
573  cout << "Impulsive start flag: "
575 
576 
577  // Set physical parameters:
578 
579  // Womersly number = Reynolds number (St = 1)
582 
583  //Horizontal length of domain
584  double lx = 1.0;
585 
586  //Vertical length of domain
587  double ly = 1.0;
588 
589  // Number of elements in x-direction
590  unsigned nx=5;
591 
592  // Number of elements in y-direction
593  unsigned ny=10;
594 
595  // Solve with Crouzeix-Raviart elements
596  {
597  // Set up doc info
598  DocInfo doc_info;
599  doc_info.number()=0;
600  doc_info.set_directory("RESLT_CR");
601 
602  //Set up problem
604  problem(nx,ny,lx,ly);
605 
606  // Run the unsteady simulation
607  problem.unsteady_run(doc_info);
608  }
609 
610 
611  // Solve with Taylor-Hood elements
612  {
613  // Set up doc info
614  DocInfo doc_info;
615  doc_info.number()=0;
616  doc_info.set_directory("RESLT_TH");
617 
618  //Set up problem
620  problem(nx,ny,lx,ly);
621 
622  // Run the unsteady simulation
623  problem.unsteady_run(doc_info);
624  }
625 
626 } // end of main
Rayleigh-type problem: 2D channel flow driven by traction bc.
RectangularQuadMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
void create_traction_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &surface_mesh_pt)
Create traction elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to the Mes...
unsigned Impulsive_start_flag
Flag for impulsive start: Default = start from exact time-periodic solution.
void doc_solution(DocInfo &doc_info)
Doc the solution.
int main(int argc, char *argv[])
Driver code for Rayleigh-type problem.
void actions_after_newton_solve()
Update after solve is empty.
RayleighTractionProblem(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly)
Problem constructor.
void actions_before_newton_solve()
Update before solve is empty.
unsigned Long_run_flag
Flag for long/short run: Default = perform long run.
void prescribed_traction(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Traction required at the top boundary.
void unsteady_run(DocInfo &doc_info)
Run an unsteady simulation.
void get_exact_u(const double &t, const Vector< double > &x, Vector< double > &u)
Exact solution of the problem as a vector.
double ReSt
Womersley = Reynolds times Strouhal.
void set_initial_condition()
Set initial condition (incl previous timesteps) according to specified function.
void actions_before_implicit_timestep()
Actions before timestep (empty)
double Re
Reynolds number.
Mesh * Surface_mesh_pt
Pointer to the "surface" mesh.