osc_ring_macro.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 flow, driven by oscillating ring
31 // with pseudo-elasticity: The mean radius of ring is adjusted to
32 // allow conservation of volume (area).
33 
34 // Oomph-lib includes
35 #include "generic.h"
36 #include "navier_stokes.h"
37 
38 //Need to instantiate templated mesh
39 #include "meshes/quarter_circle_sector_mesh.h"
40 
41 //Include namespace containing Sarah's asymptotics
43 
44 
45 using namespace std;
46 
47 using namespace oomph;
48 
49 using namespace MathematicalConstants;
50 
51 
52 
53 ////////////////////////////////////////////////////////////////////////
54 ////////////////////////////////////////////////////////////////////////
55 ////////////////////////////////////////////////////////////////////////
56 
57 
58 //==================================================
59 /// Namespace for physical parameters
60 //==================================================
61 namespace Global_Physical_Variables
62 {
63 
64  /// Reynolds number
65  double Re=100.0; // ADJUST_PRIORITY
66 
67  /// Reynolds x Strouhal number
68  double ReSt=100.0; // ADJUST_PRIORITY
69 
70 }
71 
72 
73 
74 
75 ////////////////////////////////////////////////////////////////////////
76 ////////////////////////////////////////////////////////////////////////
77 ////////////////////////////////////////////////////////////////////////
78 
79 
80 //====================================================================
81 /// Driver for oscillating ring problem: Wall performs oscillations
82 /// that resemble eigenmodes of freely oscillating ring and drives
83 /// viscous fluid flow. Mean radius of wall is adjustable and
84 /// responds to a pressure value in the fluid to allow for
85 /// mass conservation.
86 //====================================================================
87 template<class ELEMENT>
88 class OscRingNStProblem : public Problem
89 {
90 
91 public:
92 
93  /// \short Constructor: Pass timestep and function pointer to the
94  /// solution that provides the initial conditions for the fluid
95  OscRingNStProblem(const double& dt,
96  FiniteElement::UnsteadyExactSolutionFctPt IC_fct_pt);
97 
98  /// Destructor (empty)
100 
101  /// Get pointer to wall as geometric object
102  GeomObject* wall_pt()
103  {
104  return Wall_pt;
105  }
106 
107  /// Update after solve (empty)
109 
110  /// \short Update the problem specs before solve (empty)
112 
113  /// \short Update the problem specs before checking Newton
114  /// convergence: Update the fluid mesh and re-set velocity
115  /// boundary conditions -- no slip velocity on the wall means
116  /// that the velocity on the wall is enslaved.
118  {
119  // Update the fluid mesh -- auxiliary update function for algebraic
120  // nodes automatically updates no slip condition.
121  fluid_mesh_pt()->node_update();
122  }
123 
124 
125  /// \short Update the problem specs after adaptation:
126  /// Set auxiliary update function that applies no slip on all
127  /// boundary nodes and choose fluid pressure dof that drives
128  /// the wall deformation
130  {
131  // Ring boundary: No slip; this also implies that the velocity needs
132  // to be updated in response to wall motion. This needs to be reset
133  // every time the mesh is changed -- there's no mechanism by which
134  // auxiliary update functions are copied to newly created nodes.
135  // (that's because unlike boundary conditions, they don't
136  // occur exclusively at boundaries)
137  unsigned ibound=1;
138  {
139  unsigned num_nod= fluid_mesh_pt()->nboundary_node(ibound);
140  for (unsigned inod=0;inod<num_nod;inod++)
141  {
142  fluid_mesh_pt()->boundary_node_pt(ibound,inod)->
143  set_auxiliary_node_update_fct_pt(
144  FSI_functions::apply_no_slip_on_moving_wall);
145  }
146  }
147 
148  // Set the reference pressure as the constant pressure in element 0
149  dynamic_cast<PseudoBucklingRingElement*>(Wall_pt)
150  ->set_reference_pressure_pt(fluid_mesh_pt()->element_pt(0)
151  ->internal_data_pt(0));
152  }
153 
154  /// Run the time integration for ntsteps steps
155  void unsteady_run(const unsigned &ntsteps, const bool& restarted,
156  DocInfo& doc_info);
157 
158  /// \short Set initial condition (incl previous timesteps) according
159  /// to specified function.
160  void set_initial_condition();
161 
162  /// Doc the solution
163  void doc_solution(DocInfo& doc_info);
164 
165  /// Access function for the fluid mesh
166  MacroElementNodeUpdateRefineableQuarterCircleSectorMesh<ELEMENT>* fluid_mesh_pt()
167  {
168  return Fluid_mesh_pt;
169  }
170 
171  /// \short Dump problem data.
172  void dump_it(ofstream& dump_file, DocInfo doc_info);
173 
174  /// \short Read problem data.
175  void restart(ifstream& restart_file);
176 
177 private:
178 
179  /// Write header for trace file
180  void write_trace_file_header();
181 
182  /// Function pointer to set the intial condition
183  FiniteElement::UnsteadyExactSolutionFctPt IC_Fct_pt;
184 
185  /// Pointer to wall
186  GeomObject* Wall_pt;
187 
188  /// Pointer to fluid mesh
189  MacroElementNodeUpdateRefineableQuarterCircleSectorMesh<ELEMENT>* Fluid_mesh_pt;
190 
191  /// Pointer to wall mesh (contains only a single GeneralisedElement)
192  Mesh* Wall_mesh_pt;
193 
194  /// Trace file
195  ofstream Trace_file;
196 
197  /// Pointer to node on coarsest mesh on which velocity is traced
198  Node* Veloc_trace_node_pt;
199 
200  /// \short Pointer to node in symmetry plane on coarsest mesh at
201  /// which velocity is traced
202  Node* Sarah_veloc_trace_node_pt;
203 
204 };
205 
206 
207 //========================================================================
208 /// Constructor: Pass (constant) timestep and function pointer to the solution
209 /// that provides the initial conditions for the fluid.
210 //========================================================================
211 template<class ELEMENT>
213 FiniteElement::UnsteadyExactSolutionFctPt IC_fct_pt) : IC_Fct_pt(IC_fct_pt)
214 {
215 
216  // Period of oscillation
217  double T=1.0;
218 
219  //Allocate the timestepper
220  add_time_stepper_pt(new BDF<4>);
221 
222  // Initialise timestep -- also sets the weights for all timesteppers
223  // in the problem.
224  initialise_dt(dt);
225 
226  // Parameters for pseudo-buckling ring
227  double eps_buckl=0.1; // ADJUST_PRIORITY
228  double ampl_ratio=-0.5; // ADJUST_PRIORITY
229  unsigned n_buckl=2; // ADJUST_PRIORITY
230  double r_0=1.0;
231 
232  // Build wall geometric element
233  Wall_pt=new PseudoBucklingRingElement(eps_buckl,ampl_ratio,n_buckl,r_0,T,
234  time_stepper_pt());
235 
236  // Fluid mesh is suspended from wall between these two Lagrangian
237  // coordinates:
238  double xi_lo=0.0;
239  double xi_hi=2.0*atan(1.0);
240 
241  // Fractional position of dividing line for two outer blocks in fluid mesh
242  double fract_mid=0.5;
243 
244  // Build fluid mesh
245  Fluid_mesh_pt=new MacroElementNodeUpdateRefineableQuarterCircleSectorMesh<ELEMENT>(
246  Wall_pt,xi_lo,fract_mid,xi_hi,time_stepper_pt());
247 
248  // Set error estimator
249  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
250  Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
251 
252  // Fluid mesh is first sub-mesh
253  add_sub_mesh(Fluid_mesh_pt);
254 
255  // Build wall mesh
256  Wall_mesh_pt=new Mesh;
257 
258  // Wall mesh is completely empty: Add Wall element in its GeneralisedElement
259  // incarnation
260  Wall_mesh_pt->add_element_pt(dynamic_cast<GeneralisedElement*>(Wall_pt));
261 
262  // Wall mesh is second sub-mesh
263  add_sub_mesh(Wall_mesh_pt);
264 
265  // Combine all submeshes into a single Mesh
266  build_global_mesh();
267 
268  // Extract pointer to node at center of mesh (this node exists
269  // at all refinement levels and can be used to doc continuous timetrace
270  // of velocities)
271  unsigned nnode=fluid_mesh_pt()->finite_element_pt(0)->nnode();
272  Veloc_trace_node_pt=fluid_mesh_pt()->finite_element_pt(0)->node_pt(nnode-1);
273 
274  // Extract pointer to node in symmetry plane (this node exists
275  // at all refinement levels and can be used to doc continuous timetrace
276  // of velocities)
277  unsigned nnode_1d=dynamic_cast<ELEMENT*>(
278  fluid_mesh_pt()->finite_element_pt(0))->nnode_1d();
280  finite_element_pt(0)->node_pt(nnode_1d-1);
281 
282  // The "pseudo-elastic" wall element is "loaded" by a pressure.
283  // Use the "constant" pressure component in Crouzeix Raviart
284  // fluid element as that pressure.
285  dynamic_cast<PseudoBucklingRingElement*>(Wall_pt)
286  ->set_reference_pressure_pt(fluid_mesh_pt()->element_pt(0)
287  ->internal_data_pt(0));
288 
289 
290  // Set the boundary conditions for this problem:
291  //----------------------------------------------
292 
293  // All nodes are free by default -- just pin the ones that have
294  // Dirichlet conditions here.
295 
296  // Bottom boundary:
297  unsigned ibound=0;
298  {
299  unsigned num_nod= fluid_mesh_pt()->nboundary_node(ibound);
300  for (unsigned inod=0;inod<num_nod;inod++)
301  {
302  // Pin vertical velocity
303  {
304  fluid_mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
305  }
306  }
307  }
308 
309  // Ring boundary: No slip; this also implies that the velocity needs
310  // to be updated in response to wall motion
311  ibound=1;
312  {
313  unsigned num_nod=fluid_mesh_pt()->nboundary_node(ibound);
314  for (unsigned inod=0;inod<num_nod;inod++)
315  {
316  // Which node are we dealing with?
317  Node* node_pt=fluid_mesh_pt()->boundary_node_pt(ibound,inod);
318 
319  // Set auxiliary update function pointer to apply no-slip condition
320  // to velocities whenever nodal position is updated
321  node_pt->set_auxiliary_node_update_fct_pt(
322  FSI_functions::apply_no_slip_on_moving_wall);
323 
324  // Pin both velocities
325  for (unsigned i=0;i<2;i++)
326  {
327  node_pt->pin(i);
328  }
329  }
330  }
331 
332  // Left boundary:
333  ibound=2;
334  {
335  unsigned num_nod=fluid_mesh_pt()->nboundary_node(ibound);
336  for (unsigned inod=0;inod<num_nod;inod++)
337  {
338  // Pin horizontal velocity
339  {
340  fluid_mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
341  }
342  }
343  }
344 
345 
346  // Complete the build of all elements so they are fully functional
347  //----------------------------------------------------------------
348 
349  //Find number of elements in mesh
350  unsigned n_element = fluid_mesh_pt()->nelement();
351 
352  // Loop over the elements to set up element-specific
353  // things that cannot be handled by constructor
354  for(unsigned i=0;i<n_element;i++)
355  {
356  // Upcast from FiniteElement to the present element
357  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(fluid_mesh_pt()->element_pt(i));
358 
359  //Set the Reynolds number, etc
360  el_pt->re_pt() = &Global_Physical_Variables::Re;
361  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
362 
363  }
364 
365  //Attach the boundary conditions to the mesh
366  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
367 
368 
369  // Set parameters for Sarah's asymptotic solution
370  //-----------------------------------------------
371 
372  // Amplitude of the oscillation
373  SarahBL::epsilon=static_cast<PseudoBucklingRingElement*>(Wall_pt)->
374  eps_buckl();
375 
376  // Womersley number is the same as square root of Reynolds number
378 
379  // Amplitude ratio
380  SarahBL::A=static_cast<PseudoBucklingRingElement*>(Wall_pt)->ampl_ratio();
381 
382  // Buckling wavenumber
383  SarahBL::N=static_cast<PseudoBucklingRingElement*>(Wall_pt)->n_buckl_float();
384 
385  // Frequency of oscillation (period is one)
386  SarahBL::Omega=2.0*MathematicalConstants::Pi;
387 
388 }
389 
390 
391 
392 
393 
394 //========================================================================
395 /// Set initial condition: Assign previous and current values
396 /// of the velocity from the velocity field specified via
397 /// the function pointer.
398 ///
399 /// Values are assigned so that the velocities and accelerations
400 /// are correct for current time.
401 //========================================================================
402 template<class ELEMENT>
404 {
405 
406  // Elastic wall: We're starting from a given initial state in which
407  // the wall is undeformed. If set_initial_condition() is called again
408  // after mesh refinement for first timestep, this needs to be reset.
409  dynamic_cast<PseudoBucklingRingElement*>(Wall_pt)->set_R_0(1.0);
410 
411  // Backup time in global timestepper
412  double backed_up_time=time_pt()->time();
413 
414  // Past history for velocities needs to be established for t=time0-deltat, ...
415  // Then provide current values (at t=time0) which will also form
416  // the initial guess for first solve at t=time0+deltat
417 
418  // Vector of exact solution values (includes pressure)
419  Vector<double> soln(3);
420  Vector<double> x(2);
421 
422  //Find number of nodes in mesh
423  unsigned num_nod = fluid_mesh_pt()->nnode();
424 
425  // Get continuous times at previous timesteps
426  int nprev_steps=time_stepper_pt()->nprev_values();
427  Vector<double> prev_time(nprev_steps+1);
428  for (int itime=nprev_steps;itime>=0;itime--)
429  {
430  prev_time[itime]=time_pt()->time(unsigned(itime));
431  }
432 
433  // Loop over current & previous timesteps (in outer loop because
434  // the mesh might also move)
435  for (int itime=nprev_steps;itime>=0;itime--)
436  {
437  double time=prev_time[itime];
438 
439  // Set global time (because this is how the geometric object refers
440  // to continous time
441  time_pt()->time()=time;
442 
443  cout << "setting IC at time =" << time << std::endl;
444 
445  // Update the fluid mesh for this value of the continuous time
446  // (The wall object reads the continous time from the same
447  // global time object)
448  fluid_mesh_pt()->node_update();
449 
450  // Loop over the nodes to set initial guess everywhere
451  for (unsigned jnod=0;jnod<num_nod;jnod++)
452  {
453 
454  // Get nodal coordinates
455  x[0]=fluid_mesh_pt()->node_pt(jnod)->x(0);
456  x[1]=fluid_mesh_pt()->node_pt(jnod)->x(1);
457 
458  // Get intial solution
459  (*IC_Fct_pt)(time,x,soln);
460 
461  // Loop over velocity directions (velocities are in soln[0] and soln[1]).
462  for (unsigned i=0;i<2;i++)
463  {
464  fluid_mesh_pt()->node_pt(jnod)->set_value(itime,i,soln[i]);
465  }
466 
467  // Loop over coordinate directions
468  for (unsigned i=0;i<2;i++)
469  {
470  fluid_mesh_pt()->node_pt(jnod)->x(itime,i)=x[i];
471  }
472  }
473  }
474 
475  // Reset backed up time for global timestepper
476  time_pt()->time()=backed_up_time;
477 
478 }
479 
480 
481 
482 
483 
484 //========================================================================
485 /// Doc the solution
486 ///
487 //========================================================================
488 template<class ELEMENT>
489 void OscRingNStProblem<ELEMENT>::doc_solution(DocInfo& doc_info)
490 {
491 
492  cout << "Doc-ing step " << doc_info.number()
493  << " for time " << time_stepper_pt()->time_pt()->time() << std::endl;
494 
495  ofstream some_file;
496  char filename[100];
497 
498  // Number of plot points
499  unsigned npts;
500  npts=5;
501 
502 
503  // Output solution on fluid mesh
504  //-------------------------------
505  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
506  doc_info.number());
507  //some_file.precision(20);
508  some_file.open(filename);
509  unsigned nelem=fluid_mesh_pt()->nelement();
510  for (unsigned ielem=0;ielem<nelem;ielem++)
511  {
512  dynamic_cast<ELEMENT*>(fluid_mesh_pt()->element_pt(ielem))->
513  full_output(some_file,npts);
514  }
515  some_file.close();
516 
517 
518  // Plot wall posn
519  //---------------
520  sprintf(filename,"%s/Wall%i.dat",doc_info.directory().c_str(),
521  doc_info.number());
522  some_file.open(filename);
523 
524  unsigned nplot=100;
525  Vector<double > xi_wall(1);
526  Vector<double > r_wall(2);
527  for (unsigned iplot=0;iplot<nplot;iplot++)
528  {
529  xi_wall[0]=0.5*Pi*double(iplot)/double(nplot-1);
530  Wall_pt->position(xi_wall,r_wall);
531  some_file << r_wall[0] << " " << r_wall[1] << std::endl;
532  }
533  some_file.close();
534 
535 
536  // Doc Sarah's asymptotic solution
537  //--------------------------------
538  sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
539  doc_info.number());
540  some_file.open(filename);
541  fluid_mesh_pt()->output_fct(some_file,npts,
542  time_stepper_pt()->time_pt()->time(),
544  some_file.close();
545 
546 
547 
548 
549  // Get position of control point
550  //------------------------------
551  Vector<double> r(2);
552  Vector<double> xi(1);
553  xi[0]=MathematicalConstants::Pi/2.0;
554  wall_pt()->position(xi,r);
555 
556 
557 
558  // Get total volume (area) of computational domain, energies and average
559  //----------------------------------------------------------------------
560  // pressure
561  //---------
562  double area=0.0;
563  double press_int=0.0;
564  double diss=0.0;
565  double kin_en=0.0;
566  nelem=fluid_mesh_pt()->nelement();
567 
568  for (unsigned ielem=0;ielem<nelem;ielem++)
569  {
570  area+=fluid_mesh_pt()->finite_element_pt(ielem)->size();
571  press_int+=dynamic_cast<ELEMENT*>(fluid_mesh_pt()->element_pt(ielem))
572  ->pressure_integral();
573  diss+=dynamic_cast<ELEMENT*>(fluid_mesh_pt()->element_pt(ielem))->
574  dissipation();
575  kin_en+=dynamic_cast<ELEMENT*>(fluid_mesh_pt()->element_pt(ielem))->
576  kin_energy();
577  }
578 
579  // Total kinetic energy in entire domain
580  double global_kin=4.0*kin_en;
581 
582  // Max/min refinement level
583  unsigned min_level;
584  unsigned max_level;
585  fluid_mesh_pt()->get_refinement_levels(min_level,max_level);
586 
587 
588  // Total dissipation for quarter domain
589  double time=time_pt()->time();
590  double diss_asympt=SarahBL::Total_Diss_sarah(time)/4.0;
591 
592  // Asymptotic predictions for velocities at control point
593  Vector<double> x_sarah(2);
594  Vector<double> soln_sarah(3);
595  x_sarah[0]=Sarah_veloc_trace_node_pt->x(0);
596  x_sarah[1]=Sarah_veloc_trace_node_pt->x(1);
597  SarahBL::exact_soln(time,x_sarah,soln_sarah);
598 
599  // Doc
600  Trace_file << time_pt()->time()
601  << " " << r[1]
602  << " " << global_kin
603  << " " << SarahBL::Kin_energy_sarah(time_pt()->time())
604  << " " << static_cast<PseudoBucklingRingElement*>(Wall_pt)->r_0()
605  << " " << area
606  << " " << press_int/area
607  << " " << diss
608  << " " << diss_asympt
609  << " " << Veloc_trace_node_pt->x(0)
610  << " " << Veloc_trace_node_pt->x(1)
611  << " " << Veloc_trace_node_pt->value(0)
612  << " " << Veloc_trace_node_pt->value(1)
613  << " " << fluid_mesh_pt()->nelement()
614  << " " << ndof()
615  << " " << min_level
616  << " " << max_level
617  << " " << fluid_mesh_pt()->nrefinement_overruled()
618  << " " << fluid_mesh_pt()->max_error()
619  << " " << fluid_mesh_pt()->min_error()
620  << " " << fluid_mesh_pt()->max_permitted_error()
621  << " " << fluid_mesh_pt()->min_permitted_error()
622  << " " << fluid_mesh_pt()->max_keep_unrefined()
623  << " " << doc_info.number()
624  << " " << Sarah_veloc_trace_node_pt->x(0)
625  << " " << Sarah_veloc_trace_node_pt->x(1)
626  << " " << Sarah_veloc_trace_node_pt->value(0)
627  << " " << Sarah_veloc_trace_node_pt->value(1)
628  << " " << x_sarah[0]
629  << " " << x_sarah[1]
630  << " " << soln_sarah[0]
631  << " " << soln_sarah[1]
632  << " "
633  << static_cast<PseudoBucklingRingElement*>(Wall_pt)->r_0()-1.0
634  << std::endl;
635 
636 
637  // Output fluid solution on coarse-ish mesh
638  //-----------------------------------------
639 
640  // Extract all elements from quadtree representation
641  Vector<Tree*> all_element_pt;
642  fluid_mesh_pt()->forest_pt()->
643  stick_all_tree_nodes_into_vector(all_element_pt);
644 
645  // Build a coarse mesh
646  Mesh* coarse_mesh_pt = new Mesh();
647 
648  //Loop over all elements and check if their refinement level is
649  //equal to the mesh's minimum refinement level
650  nelem=all_element_pt.size();
651  for (unsigned ielem=0;ielem<nelem;ielem++)
652  {
653  Tree* el_pt=all_element_pt[ielem];
654  if (el_pt->level()==min_level)
655  {
656  coarse_mesh_pt->add_element_pt(el_pt->object_pt());
657  }
658  }
659 
660  // Output fluid solution on coarse mesh
661  sprintf(filename,"%s/coarse_soln%i.dat",doc_info.directory().c_str(),
662  doc_info.number());
663  some_file.open(filename);
664  nelem=coarse_mesh_pt->nelement();
665  for (unsigned ielem=0;ielem<nelem;ielem++)
666  {
667  dynamic_cast<ELEMENT*>(coarse_mesh_pt->element_pt(ielem))->
668  full_output(some_file,npts);
669  }
670  some_file.close();
671 
672  // Write restart file
673  sprintf(filename,"%s/restart%i.dat",doc_info.directory().c_str(),
674  doc_info.number());
675  some_file.open(filename);
676  dump_it(some_file,doc_info);
677  some_file.close();
678 
679 }
680 
681 
682 
683 
684 
685 
686 //========================================================================
687 /// Dump the solution
688 //========================================================================
689 template<class ELEMENT>
690 void OscRingNStProblem<ELEMENT>::dump_it(ofstream& dump_file,DocInfo doc_info)
691 {
692 
693  // Dump refinement status of mesh
694  //fluid_mesh_pt()->dump_refinement(dump_file);
695 
696  // Call generic dump()
697  Problem::dump(dump_file);
698 
699 }
700 
701 
702 
703 //========================================================================
704 /// Read solution from disk
705 //========================================================================
706 template<class ELEMENT>
707 void OscRingNStProblem<ELEMENT>::restart(ifstream& restart_file)
708 {
709  // Refine fluid mesh according to the instructions in restart file
710  //fluid_mesh_pt()->refine(restart_file);
711 
712  // Rebuild the global mesh
713  //rebuild_global_mesh();
714 
715  // Read generic problem data
716  Problem::read(restart_file);
717 
718 // // Complete build of all elements so they are fully functional
719 // finish_problem_setup();
720 
721  //Assign equation numbers
722  //cout <<"\nNumber of equations: " << assign_eqn_numbers()
723  // << std::endl<< std::endl;
724 }
725 
726 
727 
728 //========================================================================
729 /// Driver for timestepping the problem: Fixed timestep but
730 /// guaranteed spatial accuracy. Beautiful, innit?
731 ///
732 //========================================================================
733 template<class ELEMENT>
734 void OscRingNStProblem<ELEMENT>::unsteady_run(const unsigned& ntsteps,
735  const bool& restarted,
736  DocInfo& doc_info)
737 {
738 
739  // Open trace file
740  char filename[100];
741  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
742  Trace_file.open(filename);
743 
744  // Max. number of adaptations per solve
745  unsigned max_adapt;
746 
747  // Max. number of adaptations per solve
748  if (restarted)
749  {
750  max_adapt=0;
751  }
752  else
753  {
754  max_adapt=1;
755  }
756 
757  // Max. and min. error for adaptive refinement/unrefinement
758  fluid_mesh_pt()->max_permitted_error()= 0.5e-2;
759  fluid_mesh_pt()->min_permitted_error()= 0.5e-3;
760 
761  // Don't allow refinement to drop under given level
762  fluid_mesh_pt()->min_refinement_level()=2;
763 
764  // Don't allow refinement beyond given level
765  fluid_mesh_pt()->max_refinement_level()=6;
766 
767  // Don't bother adapting the mesh if no refinement is required
768  // and if less than ... elements are to be merged.
769  fluid_mesh_pt()->max_keep_unrefined()=20;
770 
771 
772  // Get max/min refinement levels in mesh
773  unsigned min_refinement_level;
774  unsigned max_refinement_level;
775  fluid_mesh_pt()->get_refinement_levels(min_refinement_level,
776  max_refinement_level);
777 
778  cout << "\n Initial mesh: min/max refinement levels: "
779  << min_refinement_level << " " << max_refinement_level << std::endl << std::endl;
780 
781  // Doc refinement targets
782  fluid_mesh_pt()->doc_adaptivity_targets(cout);
783 
784  // Write header for trace file
785  write_trace_file_header();
786 
787  // Doc initial condition
788  doc_solution(doc_info);
789  doc_info.number()++;
790 
791  // Switch off doc during solve
792  doc_info.disable_doc();
793 
794  // If we set first to true, then initial guess will be re-assigned
795  // after mesh adaptation. Don't want this if we've done a restart.
796  bool first;
797  bool shift;
798  if (restarted)
799  {
800  first=false;
801  shift=false;
802  // Move time back by dt to make sure we're re-solving the read-in solution
803  time_pt()->time()-=time_pt()->dt();
804  }
805  else
806  {
807  first=true;
808  shift=true;
809  }
810 
811  //Take the first fixed timestep with specified tolerance for Newton solver
812  double dt=time_pt()->dt();
813  unsteady_newton_solve(dt,max_adapt,first,shift);
814 
815  // Switch doc back on
816  doc_info.enable_doc();
817 
818  // Doc initial solution
819  doc_solution(doc_info);
820  doc_info.number()++;
821 
822  // Now do normal run; allow for one mesh adaptation per timestep
823  max_adapt=1;
824 
825  //Loop over the remaining timesteps
826  for(unsigned t=2;t<=ntsteps;t++)
827  {
828 
829  // Switch off doc during solve
830  doc_info.disable_doc();
831 
832  //Take fixed timestep
833  first=false;
834  unsteady_newton_solve(dt,max_adapt,first);
835 
836  // Switch doc back on
837  doc_info.enable_doc();
838 
839  // Doc solution
840  //if (icount%10==0)
841  {
842  doc_solution(doc_info);
843  doc_info.number()++;
844  }
845  }
846 
847  // Close trace file
848  Trace_file.close();
849 
850 }
851 
852 
853 
854 
855 //========================================================================
856 /// Write trace file header
857 //========================================================================
858 template<class ELEMENT>
860 {
861 
862  // Doc parameters in trace file
863  Trace_file << "# err_max " << fluid_mesh_pt()->max_permitted_error() << std::endl;
864  Trace_file << "# err_min " << fluid_mesh_pt()->min_permitted_error() << std::endl;
865  Trace_file << "# Re " << Global_Physical_Variables::Re << std::endl;
866  Trace_file << "# St " << Global_Physical_Variables::ReSt/
867  Global_Physical_Variables::Re << std::endl;
868  Trace_file << "# dt " << time_stepper_pt()->time_pt()->dt() << std::endl;
869  Trace_file << "# initial # elements " << mesh_pt()->nelement() << std::endl;
870  Trace_file << "# min_refinement_level "
871  << fluid_mesh_pt()->min_refinement_level() << std::endl;
872  Trace_file << "# max_refinement_level "
873  << fluid_mesh_pt()->max_refinement_level() << std::endl;
874 
875 
876 
877  Trace_file << "VARIABLES=\"time\",\"V_c_t_r_l\",\"e_k_i_n\"";
878  Trace_file << ",\"e_k_i_n_(_a_s_y_m_p_t_)\",\"R_0\",\"Area\"" ;
879  Trace_file << ",\"Average pressure\",\"Total dissipation (quarter domain)\"";
880  Trace_file << ",\"Asymptotic dissipation (quarter domain)\"";
881  Trace_file << ",\"x<sub>1</sub><sup>(track)</sup>\"";
882  Trace_file << ",\"x<sub>2</sub><sup>(track)</sup>\"";
883  Trace_file << ",\"u<sub>1</sub><sup>(track)</sup>\"";
884  Trace_file << ",\"u<sub>2</sub><sup>(track)</sup>\"";
885  Trace_file << ",\"N<sub>element</sub>\"";
886  Trace_file << ",\"N<sub>dof</sub>\"";
887  Trace_file << ",\"max. refinement level\"";
888  Trace_file << ",\"min. refinement level\"";
889  Trace_file << ",\"# of elements whose refinement was over-ruled\"";
890  Trace_file << ",\"max. error\"";
891  Trace_file << ",\"min. error\"";
892  Trace_file << ",\"max. permitted error\"";
893  Trace_file << ",\"min. permitted error\"";
894  Trace_file << ",\"max. permitted # of unrefined elements\"";
895  Trace_file << ",\"doc number\"";
896  Trace_file << ",\"x<sub>1</sub><sup>(track2 FE)</sup>\"";
897  Trace_file << ",\"x<sub>2</sub><sup>(track2 FE)</sup>\"";
898  Trace_file << ",\"u<sub>1</sub><sup>(track2 FE)</sup>\"";
899  Trace_file << ",\"u<sub>2</sub><sup>(track2 FE)</sup>\"";
900  Trace_file << ",\"x<sub>1</sub><sup>(track2 Sarah)</sup>\"";
901  Trace_file << ",\"x<sub>2</sub><sup>(track2 Sarah)</sup>\"";
902  Trace_file << ",\"u<sub>1</sub><sup>(track2 Sarah)</sup>\"";
903  Trace_file << ",\"u<sub>2</sub><sup>(track2 Sarah)</sup>\"";
904  Trace_file << ",\"R<sub>0</sub>-1\"";
905  Trace_file << std::endl;
906 
907 }
908 
909 
910 
911 
912 ////////////////////////////////////////////////////////////////////////
913 ////////////////////////////////////////////////////////////////////////
914 ////////////////////////////////////////////////////////////////////////
915 
916 
917 
918 
919 
920 
921 
922 //========================================================================
923 /// Demonstrate how to solve OscRingNSt problem in deformable domain
924 /// with mesh adaptation
925 //========================================================================
926 int main(int argc, char* argv[])
927 {
928  // Store command line arguments
929  CommandLineArgs::setup(argc,argv);
930 
931  //Do a certain number of timesteps per period
932  unsigned nstep_per_period=40; // 80; // ADJUST_PRIORITY
933  unsigned nperiod=3;
934 
935  // Work out total number of steps and timestep
936  unsigned nstep=nstep_per_period*nperiod;
937  double dt=1.0/double(nstep_per_period);
938 
939  // Set up the problem: Pass timestep and Sarah's asymptotic solution for
940  // generation of initial condition
942  problem(dt,&SarahBL::full_exact_soln);
943 
944 
945  // Restart?
946  //---------
947  bool restarted=false;
948 
949  // Pointer to restart file
950  ifstream* restart_file_pt=0;
951 
952  // No restart
953  //-----------
954  if (CommandLineArgs::Argc!=2)
955  {
956  cout << "No restart" << std::endl;
957  restarted=false;
958 
959  // Refine uniformly
960  problem.refine_uniformly();
961  problem.refine_uniformly();
962 
963  // Set initial condition on uniformly refined domain (if this isn't done
964  // then the mesh contains the interpolation of the initial guess
965  // on the original coarse mesh -- if the nodal values were zero in
966  // the interior and nonzero on the boundary, then the the waviness of
967  // of the interpolated i.g. between the nodes on the coarse mesh
968  // gets transferred onto the fine mesh where we can do better
969  problem.set_initial_condition();
970 
971  }
972 
973  // Restart
974  //--------
975  else if (CommandLineArgs::Argc==2)
976  {
977  restarted=true;
978 
979  // Open restart file
980  restart_file_pt=new ifstream(CommandLineArgs::Argv[1],ios_base::in);
981  if (restart_file_pt!=0)
982  {
983  cout << "Have opened " << CommandLineArgs::Argv[1] <<
984  " for restart. " << std::endl;
985  }
986  else
987  {
988  std::ostringstream error_stream;
989  error_stream << "ERROR while trying to open "
990  << CommandLineArgs::Argv[2]
991  << " for restart." << std::endl;
992 
993  throw OomphLibError(error_stream.str(),
994  OOMPH_CURRENT_FUNCTION,
995  OOMPH_EXCEPTION_LOCATION);
996  }
997  // Do the actual restart
998  problem.restart(*restart_file_pt);
999  }
1000 
1001 
1002  // Two command line arguments: do validation run
1003  if (CommandLineArgs::Argc==3)
1004  {
1005  nstep=3;
1006  cout << "Only doing nstep steps for validation: " << nstep << std::endl;
1007  }
1008 
1009 
1010 
1011 
1012  // Setup labels for output
1013  DocInfo doc_info;
1014 
1015  // Output directory
1016  doc_info.set_directory("RESLT");
1017 
1018 
1019  // Do unsteady run of the problem for nstep steps
1020  //-----------------------------------------------
1021  problem.unsteady_run(nstep,restarted,doc_info);
1022 
1023 
1024  // Validate the restart procedure
1025  //-------------------------------
1026  if (CommandLineArgs::Argc==3)
1027  {
1028 
1029  // Build problem and restart
1030 
1031  // Set up the problem: Pass timestep and Sarah's asymptotic solution for
1032  // generation of initial condition
1033  OscRingNStProblem<MacroElementNodeUpdateElement<
1034  RefineableQCrouzeixRaviartElement<2> > >
1035  restarted_problem(dt,&SarahBL::full_exact_soln);
1036 
1037  // Setup labels for output
1038  DocInfo restarted_doc_info;
1039 
1040  // Output directory
1041  restarted_doc_info.set_directory("RESLT_restarted");
1042 
1043  // Step number
1044  restarted_doc_info.number()=0;
1045 
1046  // Copy by performing a restart from old problem
1047  restart_file_pt=new ifstream("RESLT/restart2.dat");
1048 
1049  // Do the actual restart
1050  restarted_problem.restart(*restart_file_pt);
1051 
1052  // Do unsteady run of the problem for one step
1053  unsigned nstep=2;
1054  bool restarted=true;
1055  restarted_problem.unsteady_run(nstep,restarted,restarted_doc_info);
1056 
1057  }
1058 
1059 }
1060 
1061 
1062 
1063 
1064 
1065 
1066 
1067 
1068 
1069 
OscRingNStProblem(const double &dt, FiniteElement::UnsteadyExactSolutionFctPt IC_fct_pt)
Constructor: Pass timestep and function pointer to the solution that provides the initial conditions ...
void actions_after_newton_solve()
Update after solve (empty)
void doc_solution(DocInfo &doc_info)
Doc the solution.
AlgebraicRefineableQuarterCircleSectorMesh< ELEMENT > * fluid_mesh_pt()
Access function for the fluid mesh.
Node * Sarah_veloc_trace_node_pt
Pointer to node in symmetry plane on coarsest mesh at which velocity is traced.
void unsteady_run(const unsigned &ntsteps, const bool &restarted, DocInfo &doc_info)
Run the time integration for ntsteps steps.
void dump_it(ofstream &dump_file, DocInfo doc_info)
Dump problem data.
double ReSt
Reynolds x Strouhal number.
Definition: fsi_osc_ring.cc:81
void set_initial_condition()
Set initial condition (incl previous timesteps) according to specified function.
void actions_before_newton_solve()
Update the problem specs before solve (empty)
int main(int argc, char *argv[])
double Re
Reynolds number.
Definition: fsi_osc_ring.cc:78
void full_exact_soln(const double &time, const Vector< double > &x, Vector< double > &soln)
Full exact solution: x,y,u,v,p,du/dt,dv/dt,diss.
void restart(ifstream &restart_file)
Read problem data.
void actions_after_adapt()
Update the problem specs after adaptation: Set auxiliary update function that applies no slip on all ...
GeomObject * Wall_pt
Pointer to wall.
~OscRingNStProblem()
Destructor (empty)
MacroElementNodeUpdateRefineableQuarterCircleSectorMesh< ELEMENT > * fluid_mesh_pt()
Access function for the fluid mesh.
AlgebraicRefineableQuarterCircleSectorMesh< ELEMENT > * Fluid_mesh_pt
Pointer to fluid mesh.
MacroElementNodeUpdateRefineableQuarterCircleSectorMesh< ELEMENT > * Fluid_mesh_pt
Pointer to fluid mesh.
double Total_Diss_sarah(double t)
GeomObject * wall_pt()
Get pointer to wall as geometric object.
void write_trace_file_header()
Write header for trace file.
Mesh * Wall_mesh_pt
Pointer to wall mesh (contains only a single GeneralisedElement)
Node * Veloc_trace_node_pt
Pointer to node on coarsest mesh on which velocity is traced.
void actions_before_newton_convergence_check()
Update the problem specs before checking Newton convergence: Update the fluid mesh and re-set velocit...
void exact_soln(const double &time, const Vector< double > &x, Vector< double > &soln)
Exact solution: x,y,u,v,p.
double Kin_energy_sarah(double t)