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