osc_quarter_ellipse.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 2D Navier-Stokes problem in moving domain
31 
32 //Generic routines
33 #include "generic.h"
34 
35 // The Navier Stokes equations
36 #include "navier_stokes.h"
37 
38 // Mesh
39 #include "meshes/quarter_circle_sector_mesh.h"
40 
41 using namespace std;
42 
43 using namespace oomph;
44 
45 //============start_of_MyEllipse===========================================
46 /// \short Oscillating ellipse
47 /// \f[ x = (A + \widehat{A} \cos(2\pi t/T)) \cos(\xi) \f]
48 /// \f[ y = \frac{\sin(\xi)}{A + \widehat{A} \cos(2\pi t/T)} \f]
49 /// Note that cross-sectional area is conserved.
50 //=========================================================================
51 class MyEllipse : public GeomObject
52 {
53 
54 public:
55 
56  /// \short Constructor: Pass initial x-half axis, amplitude of x-variation,
57  /// period of oscillation and pointer to time object.
58  MyEllipse(const double& a, const double& a_hat,
59  const double& period, Time* time_pt) :
60  GeomObject(1,2), A(a), A_hat(a_hat), T(period), Time_pt(time_pt) {}
61 
62  /// Destructor: Empty
63  virtual ~MyEllipse() {}
64 
65  /// \short Current position vector to material point at
66  /// Lagrangian coordinate xi
67  void position(const Vector<double>& xi, Vector<double>& r) const
68  {
69  // Get current time:
70  double time=Time_pt->time();
71 
72  // Position vector
73  double axis=A+A_hat*cos(2.0*MathematicalConstants::Pi*time/T);
74  r[0] = axis*cos(xi[0]);
75  r[1] = (1.0/axis)*sin(xi[0]);
76  }
77 
78  /// \short Parametrised position on object: r(xi). Evaluated at
79  /// previous time level. t=0: current time; t>0: previous
80  /// time level.
81  void position(const unsigned& t, const Vector<double>& xi,
82  Vector<double>& r) const
83  {
84  // Get current time:
85  double time=Time_pt->time(t);
86 
87  // Position vector
88  double axis=A+A_hat*cos(2.0*MathematicalConstants::Pi*time/T);
89  r[0] = axis*cos(xi[0]);
90  r[1] = (1.0/axis)*sin(xi[0]);
91  }
92 
93 private:
94 
95  /// x-half axis
96  double A;
97 
98  /// Amplitude of variation in x-half axis
99  double A_hat;
100 
101  /// Period of oscillation
102  double T;
103 
104  /// Pointer to time object
105  Time* Time_pt;
106 
107 }; // end of MyEllipse
108 
109 
110 
111 ///////////////////////////////////////////////////////////////////////
112 ///////////////////////////////////////////////////////////////////////
113 ////////////////////////////////////////////////////////////////////////
114 
115 
116 
117 //===start_of_namespace=================================================
118 /// Namepspace for global parameters
119 //======================================================================
120 namespace Global_Physical_Variables
121 {
122 
123  /// Reynolds number
124  double Re=100.0;
125 
126  /// Womersley = Reynolds times Strouhal
127  double ReSt=100.0;
128 
129  /// x-Half axis length
130  double A=1.0;
131 
132  /// x-Half axis amplitude
133  double A_hat=0.1;
134 
135  /// Period of oscillations
136  double T=1.0;
137 
138  /// Exact solution of the problem as a vector containing u,v,p
139  void get_exact_u(const double& t, const Vector<double>& x, Vector<double>& u)
140  {
141  using namespace MathematicalConstants;
142 
143  // Strouhal number
144  double St = ReSt/Re;
145 
146  // Half axis
147  double a=A+A_hat*cos(2.0*Pi*t/T);
148  double adot=-2.0*A_hat*Pi*sin(2.0*Pi*t/T)/T;
149  u.resize(3);
150 
151  // Velocity solution
152  u[0]=adot*x[0]/a;
153  u[1]=-adot*x[1]/a;
154 
155  // Pressure solution
156  u[2]=(2.0*A_hat*Pi*Pi*Re*(x[0]*x[0]*St*cos(2.0*Pi*t/T)*A +
157  x[0]*x[0]*St*A_hat - x[0]*x[0]*A_hat +
158  x[0]*x[0]*A_hat*cos(2.0*Pi*t/T)*cos(2.0*Pi*t/T) -
159  x[1]*x[1]*St*cos(2.0*Pi*t/T)*A -
160  x[1]*x[1]*St*A_hat - x[1]*x[1]*A_hat +
161  x[1]*x[1]*A_hat*cos(2.0*Pi*t/T)*cos(2.0*Pi*t/T) ))
162  /(T*T*(A*A + 2.0*A*A_hat*cos(2.0*Pi*t/T) +
163  A_hat*A_hat*cos(2.0*Pi*t/T)*cos(2.0*Pi*t/T) ));
164  }
165 
166 
167 } // end of namespace
168 
169 
170 
171 
172 ///////////////////////////////////////////////////////////////////////
173 ///////////////////////////////////////////////////////////////////////
174 ////////////////////////////////////////////////////////////////////////
175 
176 
177 
178 //=====start_of_problem_class=========================================
179 /// Navier-Stokes problem in an oscillating ellipse domain.
180 //====================================================================
181 template<class ELEMENT, class TIMESTEPPER>
182 class OscEllipseProblem : public Problem
183 {
184 
185 public:
186 
187  /// Constructor
189 
190  /// Destructor (empty)
192 
193  /// Update the problem specs after solve (empty)
195 
196  /// \short Update problem specs before solve (empty)
198 
199  /// Actions before adapt (empty)
201 
202  /// Actions after adaptation, pin relevant pressures
204  {
205  // Unpin all pressure dofs
206  RefineableNavierStokesEquations<2>::
207  unpin_all_pressure_dofs(mesh_pt()->element_pt());
208 
209  // Pin redundant pressure dofs
210  RefineableNavierStokesEquations<2>::
211  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
212 
213  // Now set the first pressure dof in the first element to 0.0
214  fix_pressure(0,0,0.0);
215 
216  } // end of actions_after_adapt
217 
218 
219  /// \short Update the problem specs before next timestep
221  {
222  // Update the domain shape
223  mesh_pt()->node_update();
224 
225  // Ring boundary: No slip; this implies that the velocity needs
226  // to be updated in response to wall motion
227  unsigned ibound=1;
228  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
229  for (unsigned inod=0;inod<num_nod;inod++)
230  {
231  // Which node are we dealing with?
232  Node* node_pt=mesh_pt()->boundary_node_pt(ibound,inod);
233 
234  // Apply no slip
235  FSI_functions::apply_no_slip_on_moving_wall(node_pt);
236  }
237  }
238 
239  /// Update the problem specs after timestep (empty)
241 
242  /// Doc the solution
243  void doc_solution(DocInfo& doc_info);
244 
245  /// Timestepping loop
246  void unsteady_run(DocInfo& doc_info);
247 
248  /// \short Set initial condition
249  void set_initial_condition();
250 
251 private:
252 
253  ///Fix pressure in element e at pressure dof pdof and set to pvalue
254  void fix_pressure(const unsigned &e, const unsigned &pdof,
255  const double &pvalue)
256  {
257  //Cast to proper element and fix pressure
258  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
259  fix_pressure(pdof,pvalue);
260  } // end_of_fix_pressure
261 
262  /// Pointer to GeomObject that specifies the domain bondary
263  GeomObject* Wall_pt;
264 
265 }; // end of problem_class
266 
267 
268 
269 //========start_of_constructor============================================
270 /// Constructor for Navier-Stokes problem on an oscillating ellipse domain.
271 //========================================================================
272 template<class ELEMENT, class TIMESTEPPER>
274 {
275 
276  //Create the timestepper and add it to the problem
277  add_time_stepper_pt(new TIMESTEPPER);
278 
279  // Setup mesh
280  //-----------
281 
282  // Build geometric object that forms the curvilinear domain boundary:
283  // an oscillating ellipse
284 
285  // Half axes
287 
288  // Variations of half axes
290 
291  // Period of the oscillation
292  double period=Global_Physical_Variables::T;
293 
294  // Create GeomObject that specifies the domain bondary
295  Wall_pt=new MyEllipse(a,a_hat,period,Problem::time_pt());
296 
297 
298  // Start and end coordinates of curvilinear domain boundary on ellipse
299  double xi_lo=0.0;
300  double xi_hi=MathematicalConstants::Pi/2.0;
301 
302  // Now create the mesh. Separating line between the two
303  // elements next to the curvilinear boundary is located half-way
304  // along the boundary.
305  double fract_mid=0.5;
306  Problem::mesh_pt() = new RefineableQuarterCircleSectorMesh<ELEMENT>(
307  Wall_pt,xi_lo,fract_mid,xi_hi,time_stepper_pt());
308 
309  // Set error estimator NOT NEEDED IN CURRENT PROBLEM SINCE
310  // WE'RE ONLY REFINING THE MESH UNIFORMLY
311  //Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
312  //mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
313 
314 
315  // Fluid boundary conditions
316  //--------------------------
317  // Ring boundary: No slip; this also implies that the velocity needs
318  // to be updated in response to wall motion
319  unsigned ibound=1;
320  {
321  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
322  for (unsigned inod=0;inod<num_nod;inod++)
323  {
324 
325  // Pin both velocities
326  for (unsigned i=0;i<2;i++)
327  {
328  mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
329  }
330  }
331  } // end boundary 1
332 
333  // Bottom boundary:
334  ibound=0;
335  {
336  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
337  for (unsigned inod=0;inod<num_nod;inod++)
338  {
339  // Pin vertical velocity
340  {
341  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
342  }
343  }
344  } // end boundary 0
345 
346  // Left boundary:
347  ibound=2;
348  {
349  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
350  for (unsigned inod=0;inod<num_nod;inod++)
351  {
352  // Pin horizontal velocity
353  {
354  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
355  }
356  }
357  } // end boundary 2
358 
359 
360  // Complete the build of all elements so they are fully functional
361  //----------------------------------------------------------------
362 
363  // Find number of elements in mesh
364  unsigned n_element = mesh_pt()->nelement();
365 
366  // Loop over the elements to set up element-specific
367  // things that cannot be handled by constructor
368  for(unsigned i=0;i<n_element;i++)
369  {
370  // Upcast from FiniteElement to the present element
371  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
372 
373  //Set the Reynolds number, etc
374  el_pt->re_pt() = &Global_Physical_Variables::Re;
375  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
376  }
377 
378  // Pin redundant pressure dofs
379  RefineableNavierStokesEquations<2>::
380  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
381 
382  // Now set the first pressure dof in the first element to 0.0
383  fix_pressure(0,0,0.0);
384 
385  // Do equation numbering
386  cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
387 
388 } // end of constructor
389 
390 
391 //======================start_of_set_initial_condition====================
392 /// \short Set initial condition: Assign previous and current values
393 /// from exact solution.
394 //========================================================================
395 template<class ELEMENT,class TIMESTEPPER>
397 {
398  // Backup time in global timestepper
399  double backed_up_time=time_pt()->time();
400 
401  // Past history for velocities must be established for t=time0-deltat, ...
402  // Then provide current values (at t=time0) which will also form
403  // the initial guess for first solve at t=time0+deltat
404 
405  // Vector of exact solution value
406  Vector<double> soln(3);
407  Vector<double> x(2);
408 
409  //Find number of nodes in mesh
410  unsigned num_nod = mesh_pt()->nnode();
411 
412  // Get continuous times at previous timesteps
413  int nprev_steps=time_stepper_pt()->nprev_values();
414  Vector<double> prev_time(nprev_steps+1);
415  for (int itime=nprev_steps;itime>=0;itime--)
416  {
417  prev_time[itime]=time_pt()->time(unsigned(itime));
418  }
419 
420  // Loop over current & previous timesteps (in outer loop because
421  // the mesh also moves!)
422  for (int itime=nprev_steps;itime>=0;itime--)
423  {
424  double time=prev_time[itime];
425 
426  // Set global time (because this is how the geometric object refers
427  // to continous time )
428  time_pt()->time()=time;
429 
430  cout << "setting IC at time =" << time << std::endl;
431 
432  // Update the mesh for this value of the continuous time
433  // (The wall object reads the continous time from the
434  // global time object)
435  mesh_pt()->node_update();
436 
437  // Loop over the nodes to set initial guess everywhere
438  for (unsigned jnod=0;jnod<num_nod;jnod++)
439  {
440  // Get nodal coordinates
441  x[0]=mesh_pt()->node_pt(jnod)->x(0);
442  x[1]=mesh_pt()->node_pt(jnod)->x(1);
443 
444  // Get exact solution (unsteady stagnation point flow)
446 
447  // Assign solution
448  mesh_pt()->node_pt(jnod)->set_value(itime,0,soln[0]);
449  mesh_pt()->node_pt(jnod)->set_value(itime,1,soln[1]);
450 
451  // Loop over coordinate directions
452  for (unsigned i=0;i<2;i++)
453  {
454  mesh_pt()->node_pt(jnod)->x(itime,i)=x[i];
455  }
456  }
457  } // end of loop over previous timesteps
458 
459  // Reset backed up time for global timestepper
460  time_pt()->time()=backed_up_time;
461 
462 } // end of set_initial_condition
463 
464 
465 
466 //=======start_of_doc_solution============================================
467 /// Doc the solution
468 //========================================================================
469 template<class ELEMENT, class TIMESTEPPER>
471 {
472  ofstream some_file;
473  char filename[100];
474 
475  // Number of plot points
476  unsigned npts;
477  npts=5;
478 
479  // Output solution
480  //-----------------
481  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
482  doc_info.number());
483  some_file.open(filename);
484  mesh_pt()->output(some_file,npts);
485  some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
486  << time_pt()->time() << "\"";
487  some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
488  some_file << "1" << std::endl;
489  some_file << "2" << std::endl;
490  some_file << " 0 0" << std::endl;
491  some_file << time_pt()->time()*20.0 << " 0" << std::endl;
492 
493  // Write dummy zones that force tecplot to keep the axis limits constant
494  // while the domain is moving.
495  some_file << "ZONE I=2,J=2" << std::endl;
496  some_file << "0.0 0.0 -0.65 -0.65 -200.0" << std::endl;
497  some_file << "1.15 0.0 -0.65 -0.65 -200.0" << std::endl;
498  some_file << "0.0 1.15 -0.65 -0.65 -200.0" << std::endl;
499  some_file << "1.15 1.15 -0.65 -0.65 -200.0" << std::endl;
500  some_file << "ZONE I=2,J=2" << std::endl;
501  some_file << "0.0 0.0 0.65 0.65 300.0" << std::endl;
502  some_file << "1.15 0.0 0.65 0.65 300.0" << std::endl;
503  some_file << "0.0 1.15 0.65 0.65 300.0" << std::endl;
504  some_file << "1.15 1.15 0.65 0.65 300.0" << std::endl;
505 
506  some_file.close();
507 
508  // Output exact solution
509  //----------------------
510  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
511  doc_info.number());
512  some_file.open(filename);
513  mesh_pt()->output_fct(some_file,npts,time_pt()->time(),
515  some_file.close();
516 
517  // Doc error
518  //----------
519  double error,norm;
520  sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
521  doc_info.number());
522  some_file.open(filename);
523  mesh_pt()->compute_error(some_file,
525  time_pt()->time(),
526  error,norm);
527  some_file.close();
528 
529 
530  // Doc solution and error
531  //-----------------------
532  cout << "error: " << error << std::endl;
533  cout << "norm : " << norm << std::endl << std::endl;
534 
535 
536  // Plot wall posn
537  //---------------
538  sprintf(filename,"%s/Wall%i.dat",doc_info.directory().c_str(),
539  doc_info.number());
540  some_file.open(filename);
541 
542  unsigned nplot=100;
543  for (unsigned iplot=0;iplot<nplot;iplot++)
544  {
545  Vector<double> xi_wall(1), r_wall(2);
546  xi_wall[0]=0.5*MathematicalConstants::Pi*double(iplot)/double(nplot-1);
547  Wall_pt->position(xi_wall,r_wall);
548  some_file << r_wall[0] << " " << r_wall[1] << std::endl;
549  }
550  some_file.close();
551 
552  // Increment number of doc
553  doc_info.number()++;
554 
555 } // end of doc_solution
556 
557 
558 //=======start_of_unsteady_run============================================
559 /// Unsteady run
560 //========================================================================
561 template<class ELEMENT, class TIMESTEPPER>
563 {
564 
565  // Specify duration of the simulation
566  double t_max=3.0;
567 
568  // Initial timestep
569  double dt=0.025;
570 
571  // Initialise timestep
572  initialise_dt(dt);
573 
574  // Set initial conditions.
575  set_initial_condition();
576 
577  // Alternative initial conditions: impulsive start; see exercise.
578  //assign_initial_values_impulsive();
579 
580  // find number of steps
581  unsigned nstep = unsigned(t_max/dt);
582 
583  // If validation: Reduce number of timesteps performed and
584  // use coarse-ish mesh
585  if (CommandLineArgs::Argc>1)
586  {
587  nstep=2;
588  refine_uniformly();
589  cout << "validation run" << std::endl;
590  }
591  else
592  {
593  // Refine the mesh three times, to resolve the pressure distribution
594  // (the velocities could be represented accurately on a much coarser mesh).
595  refine_uniformly();
596  refine_uniformly();
597  refine_uniformly();
598  }
599 
600  // Output solution initial
601  doc_solution(doc_info);
602 
603  // Timestepping loop
604  for (unsigned istep=0;istep<nstep;istep++)
605  {
606  cout << "TIMESTEP " << istep << std::endl;
607  cout << "Time is now " << time_pt()->time() << std::endl;
608 
609  // Take timestep
610  unsteady_newton_solve(dt);
611 
612  //Output solution
613  doc_solution(doc_info);
614  }
615 
616 } // end of unsteady_run
617 
618 
619 ////////////////////////////////////////////////////////////////////////
620 ////////////////////////////////////////////////////////////////////////
621 
622 
623 //======start_of_main=================================================
624 /// Driver code for unsteady Navier-Stokes flow, driven by
625 /// oscillating ellipse. If the code is executed with command line
626 /// arguments, a validation run is performed.
627 //====================================================================
628 int main(int argc, char* argv[])
629 {
630 
631  // Store command line arguments
632  CommandLineArgs::setup(argc,argv);
633 
634 
635  // Solve with Crouzeix-Raviart elements
636  {
637  // Create DocInfo object with suitable directory name for output
638  DocInfo doc_info;
639  doc_info.set_directory("RESLT_CR");
640 
641  //Set up problem
643 
644  // Run the unsteady simulation
645  problem.unsteady_run(doc_info);
646  }
647 
648  // Solve with Taylor-Hood elements
649  {
650  // Create DocInfo object with suitable directory name for output
651  DocInfo doc_info;
652  doc_info.set_directory("RESLT_TH");
653 
654  //Set up problem
656 
657  // Run the unsteady simulation
658  problem.unsteady_run(doc_info);
659  }
660 
661 
662 }; // end of main
663 
664 
665 
OscEllipseProblem()
Constructor.
void actions_after_adapt()
Actions after adaptation, pin relevant pressures.
double T
Period of oscillations.
void unsteady_run(DocInfo &doc_info)
Timestepping loop.
void set_initial_condition()
Set initial condition.
Time * Time_pt
Pointer to time object.
void actions_before_newton_solve()
Update problem specs before solve (empty)
GeomObject * Wall_pt
Pointer to GeomObject that specifies the domain bondary.
void actions_after_implicit_timestep()
Update the problem specs after timestep (empty)
Navier-Stokes problem in an oscillating ellipse domain.
void position(const Vector< double > &xi, Vector< double > &r) const
Current position vector to material point at Lagrangian coordinate xi.
~OscEllipseProblem()
Destructor (empty)
double ReSt
Womersley = Reynolds times Strouhal.
double T
Period of oscillation.
void actions_before_implicit_timestep()
Update the problem specs before next timestep.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Parametrised position on object: r(xi). Evaluated at previous time level. t=0: current time; t>0: pre...
double Re
Reynolds number.
int main(int argc, char *argv[])
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void doc_solution(DocInfo &doc_info)
Doc the solution.
virtual ~MyEllipse()
Destructor: Empty.
Oscillating ellipse Note that cross-sectional area is conserved.
MyEllipse(const double &a, const double &a_hat, const double &period, Time *time_pt)
Constructor: Pass initial x-half axis, amplitude of x-variation, period of oscillation and pointer to...
void get_exact_u(const double &t, const Vector< double > &x, Vector< double > &u)
Exact solution of the problem as a vector containing u,v,p.
double A_hat
x-Half axis amplitude
double A
x-half axis
double A
x-Half axis length
void actions_before_adapt()
Actions before adapt (empty)
double A_hat
Amplitude of variation in x-half axis.
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.