fsi_osc_ring.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 interacting with an elastic ring
31 
32 // Oomph-lib include files
33 #include "generic.h"
34 #include "navier_stokes.h"
35 #include "beam.h"
36 
37 //Need to include templated meshes, so that all functions
38 //are instantiated for our particular element types.
39 #include "meshes/quarter_circle_sector_mesh.h"
40 #include "meshes/one_d_lagrangian_mesh.h"
41 
42 using namespace std;
43 
44 using namespace oomph;
45 
46 //===========================================================================
47 /// Namespace for physical parameters
48 //===========================================================================
49 namespace Global_Physical_Variables
50 {
51 
52  // Independent parameters:
53  //------------------------
54 
55  /// Square of Womersly number (a frequency parameter)
56  double Alpha_sq=50.0;
57 
58  /// Density ratio of the solid and the fluid
59  double Density_ratio=1.0;
60 
61  /// External Pressure
62  double Pext=0.0;
63 
64  /// Poisson ratio
65  double Nu=0.49;
66 
67  /// Nondimensional thickness of the beam
68  double H=0.05;
69 
70  /// Perturbation pressure
71  double Pcos=0.0;
72 
73 
74  // Dependent parameters:
75  //----------------------
76 
77  /// Reynolds number
78  double Re;
79 
80  /// Reynolds x Strouhal number
81  double ReSt;
82 
83  /// Timescale ratio (non-dimensation density)
84  double Lambda_sq;
85 
86  /// Stress ratio
87  double Q;
88 
89  /// \short Set the parameters that are used in the code as a function
90  /// of the Womersley number, the density ratio and H
91  void set_params()
92  {
93  cout << "\n\n======================================================" <<std::endl;
94  cout << "\nSetting parameters. \n\n";
95  cout << "Prescribed: Square of Womersley number: Alpha_sq = "
96  << Alpha_sq << std::endl;
97  cout << " Density ratio: Density_ratio = "
98  << Density_ratio << std::endl;
99  cout << " Wall thickness: H = "
100  << H << std::endl;
101  cout << " Poisson ratio: Nu = "
102  << Nu << std::endl;
103  cout << " Pressure perturbation: Pcos = "
104  << Pcos << std::endl;
105 
106 
107  Q=1.0/12.0*pow(H,3)/Alpha_sq;
108  cout << "\nDependent: Stress ratio: Q = "
109  << Q << std::endl;
110 
111  Lambda_sq=1.0/12.0*pow(H,3)*Density_ratio;
112  cout << " Timescale ratio: Lambda_sq = "
113  << Lambda_sq << std::endl;
114 
115  Re=Alpha_sq;
116  cout << " Reynolds number: Re = "
117  << Re << std::endl;
118 
119  ReSt=Re;
120  cout << " Womersley number: ReSt = "
121  << ReSt << std::endl;
122  cout << "\n======================================================\n\n"
123  <<std::endl;
124 
125  }
126 
127  /// Non-FSI load function, a constant external pressure plus
128  /// a (small) sinusoidal perturbation of wavenumber two.
129  void pcos_load(const Vector<double>& xi, const Vector<double> &x,
130  const Vector<double>& N, Vector<double>& load)
131  {
132  for(unsigned i=0;i<2;i++)
133  {load[i] = (Pext - Pcos*cos(2.0*xi[0]))*N[i];}
134  }
135 
136 }
137 
138 
139 //======================================================================
140 /// FSI Ring problem: a fluid-structure interaction problem in which
141 /// a viscous fluid bounded by an initially circular beam is set into motion
142 /// by a small sinusoidal perturbation of the beam (the domain boundary).
143 //======================================================================
144 class FSIRingProblem : public Problem
145 {
146  /// \short There are very few element types that will work for this problem.
147  /// Rather than passing the element type as a template parameter to the
148  /// problem, we choose instead to use a typedef to specify the
149  /// particular element fluid used.
150  typedef AlgebraicElement<RefineableQCrouzeixRaviartElement<2> > FLUID_ELEMENT;
151 
152  /// \short Typedef to specify the solid element used
153  typedef FSIHermiteBeamElement SOLID_ELEMENT;
154 
155 public:
156 
157  /// Constructor: Number of elements in wall mesh, amplitude of the
158  /// initial wall deformation, amplitude of pcos perturbation and its duration.
159  FSIRingProblem(const unsigned &nelement_wall,
160  const double& eps_ampl, const double& pcos_initial,
161  const double& pcos_duration);
162 
163  /// Update after solve (empty)
165 
166  /// \short Update before solve (empty)
168 
169  /// \short Update the problem specs before checking Newton
170  /// convergence
172  {
173  // Update the fluid mesh -- auxiliary update function for algebraic
174  // nodes automatically updates no slip condition.
175  Fluid_mesh_pt->node_update();
176  }
177 
178  /// \short Update the problem specs after adaptation:
180  {
181  // The functions used to update the no slip boundary conditions
182  // must be set on any new nodes that have been created during the
183  // mesh adaptation process.
184  // There is no mechanism by which auxiliary update functions
185  // are copied to newly created nodes.
186  // (because, unlike boundary conditions, they don't occur exclusively
187  // at boundaries)
188 
189  // The no-slip boundary is boundary 1 of the mesh
190  // Loop over the nodes on this boundary and reset the auxilliary
191  // node update function
192  unsigned n_node = Fluid_mesh_pt->nboundary_node(1);
193  for (unsigned n=0;n<n_node;n++)
194  {
195  Fluid_mesh_pt->boundary_node_pt(1,n)->set_auxiliary_node_update_fct_pt(
196  FSI_functions::apply_no_slip_on_moving_wall);
197  }
198 
199  // (Re-)setup fsi: Work out which fluid dofs affect wall elements
200  // the correspondance between wall dofs and fluid elements is handled
201  // during the remeshing, but the "reverse" association must be done
202  // separately.
203  // We need to set up the interaction every time because the fluid element
204  // adjacent to a given solid element's integration point may have changed
205  // We pass the boundary between the fluid and solid meshes and pointers
206  // to the meshes. The interaction boundary is boundary 1 of the 2D
207  // fluid mesh.
208  FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
209  (this,1,Fluid_mesh_pt,Wall_mesh_pt);
210  }
211 
212  /// \short Doc solution: Pass number of timestep, i (we append to tracefile
213  /// after every timestep but do a full doc only at certain intervals),
214  /// DocInfo object and tracefile
215  void doc_solution(const unsigned& i, DocInfo& doc_info, ofstream& trace_file);
216 
217  /// Do dynamic run
218  void dynamic_run();
219 
220 private:
221 
222  /// Setup initial condition for both domains
223  void set_initial_condition();
224 
225  /// Setup initial condition for wall
226  void set_wall_initial_condition();
227 
228  /// Setup initial condition for fluid
229  void set_fluid_initial_condition();
230 
231  /// \short Element used for documenting displacement
233 
234  /// Pointer to wall mesh
235  OneDLagrangianMesh<SOLID_ELEMENT> *Wall_mesh_pt;
236 
237  /// Pointer to fluid mesh
238  AlgebraicRefineableQuarterCircleSectorMesh<FLUID_ELEMENT> *Fluid_mesh_pt;
239 
240  /// Pointer to geometric object that represents the undeformed wall shape
241  GeomObject* Undef_geom_pt;
242 
243  /// Pointer to wall timestepper
244  Newmark<2>* Wall_time_stepper_pt;
245 
246  /// Pointer to fluid timestepper
248 
249  /// Pointer to node on coarsest mesh on which velocity is traced
251 
252  /// Amplitude of initial deformation
253  double Eps_ampl;
254 
255  /// Initial pcos
256  double Pcos_initial;
257 
258  /// Duration of initial pcos
260 
261 };
262 
263 
264 //===============================================================
265 /// Setup initial condition: When we're done here, all variables
266 /// represent the state at the initial time.
267 //===============================================================
269 {
270 
271  cout << "Setting wall ic" << std::endl;
272  set_wall_initial_condition();
273 
274  cout << "Setting fluid ic" << std::endl;
275  set_fluid_initial_condition();
276 
277 }
278 
279 
280 //===============================================================
281 /// Setup initial condition for fluid: Impulsive start
282 //===============================================================
284 {
285 
286  // Update fluid domain: Careful!!! This also applies the no slip conditions
287  // on all nodes on the wall! Since the wall might have moved since
288  // we created the mesh; we're therefore imposing a nonzero
289  // velocity on these nodes. Must wipe this afterwards (done
290  // by setting *all* velocities to zero) otherwise we get
291  // an impulsive start from a very bizarre initial velocity
292  // field! [Yes, it took me a while to figure this out...]
293  Fluid_mesh_pt->node_update();
294 
295  // Assign initial values for the velocities;
296  // pressures don't carry a time history and can be left alone.
297 
298  //Find number of nodes in fluid mesh
299  unsigned n_node = Fluid_mesh_pt->nnode();
300 
301  // Loop over the nodes to set initial guess everywhere
302  for(unsigned n=0;n<n_node;n++)
303  {
304  // Loop over velocity directions: Impulsive initial start from
305  // zero velocity!
306  for(unsigned i=0;i<2;i++)
307  {
308  Fluid_mesh_pt->node_pt(n)->set_value(i,0.0);
309  }
310  }
311 
312  // Do an impulsive start with the assigned velocity field
313  Fluid_mesh_pt->assign_initial_values_impulsive();
314 
315 }
316 
317 
318 //===============================================================
319 /// Setup initial condition: Impulsive start either from
320 /// deformed or undeformed wall shape.
321 //===============================================================
323 {
324 
325  // Geometric object that specifies the initial conditions:
326  // A ring that is bucked in a 2-lobed mode
327  GeomObject* ic_geom_object_pt=
328  new PseudoBucklingRing(Eps_ampl,Global_Physical_Variables::H,2,2,
329  Wall_time_stepper_pt);
330 
331  // Assign period of oscillation of the geometric object
332  static_cast<PseudoBucklingRing*>(ic_geom_object_pt)->set_T(1.0);
333 
334  //Set initial time (to deform wall into max. amplitude)
335  double time=0.25;
336 
337  // Assign initial radius of the object
338  static_cast<PseudoBucklingRing*>(ic_geom_object_pt)->set_R_0(1.00);
339 
340  // Setup object that specifies the initial conditions:
341  SolidInitialCondition* IC_pt = new SolidInitialCondition(ic_geom_object_pt);
342 
343  // Assign values of positional data of all elements on wall mesh
344  // so that the wall deforms into the shape specified by IC object.
345  SolidMesh::Solid_IC_problem.set_static_initial_condition(
346  this,Wall_mesh_pt,IC_pt,time);
347 
348 }
349 
350 
351 //========================================================================
352 /// Document solution: Pass number of timestep, i; we append to trace file
353 /// at every timestep and do a full doc only after a certain number
354 /// of steps.
355 //========================================================================
356 void FSIRingProblem::doc_solution(const unsigned& i,
357  DocInfo& doc_info, ofstream& trace_file)
358 {
359 
360  // Full doc every nskip steps
361  unsigned nskip=1; // ADJUST
362 
363  // If we at an integer multiple of nskip, full documentation.
364  if (i%nskip==0)
365  {
366  doc_info.enable_doc();
367  cout << "Full doc step " << doc_info.number()
368  << " for time " << time_stepper_pt()->time_pt()->time() << std::endl;
369  }
370  //Otherwise, just output the trace file
371  else
372  {
373  doc_info.disable_doc();
374  cout << "Only trace for time "
375  << time_stepper_pt()->time_pt()->time() << std::endl;
376  }
377 
378 
379  // If we are at a full documentation step, output the fluid solution
380  if (doc_info.is_doc_enabled())
381  {
382  //Variables used in the output file.
383  ofstream some_file; char filename[100];
384  //Construct the output filename from the doc_info number and the
385  //output directory
386  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
387  doc_info.number());
388  //Open the output file
389  some_file.open(filename);
390  ///Output the solution using 5x5 plot points
391  Fluid_mesh_pt->output(some_file,5);
392  //Close the output file
393  some_file.close();
394  }
395 
396  //Temporary vector to give the local coordinate at which to document
397  //the wall displacment
398  Vector<double> s(1,1.0);
399  // Write to the trace file:
400  trace_file << time_pt()->time()
401  //Document the displacement at the end of the the chosen element
402  << " " << Doc_displacement_elem_pt->interpolated_x(s,1)
403  << " " << Veloc_trace_node_pt->x(0)
404  << " " << Veloc_trace_node_pt->x(1)
405  << " " << Veloc_trace_node_pt->value(0)
406  << " " << Veloc_trace_node_pt->value(1)
407  << " " << Fluid_mesh_pt->nelement()
408  << " " << ndof()
409  << " " << Fluid_mesh_pt->nrefinement_overruled()
410  << " " << Fluid_mesh_pt->max_error()
411  << " " << Fluid_mesh_pt->min_error()
412  << " " << Fluid_mesh_pt->max_permitted_error()
413  << " " << Fluid_mesh_pt->min_permitted_error()
414  << " " << Fluid_mesh_pt->max_keep_unrefined();
415 
416  // Output the number of the corresponding full documentation
417  // file number (or -1 if no full doc was made)
418  if (doc_info.is_doc_enabled())
419  {trace_file << " " <<doc_info.number() << " ";}
420  else {trace_file << " " <<-1 << " ";}
421 
422  //End the trace file
423  trace_file << std::endl;
424 
425  // Increment counter for full doc
426  if (doc_info.is_doc_enabled()) {doc_info.number()++;}
427 }
428 
429 //======================================================================
430 /// Constructor for FSI ring problem. Pass number of wall elements
431 /// and length of wall (in Lagrangian coordinates) amplitude of
432 /// initial deformation, pcos perturbation and duration.
433 //======================================================================
435  const double& eps_ampl, const double& pcos_initial,
436  const double& pcos_duration) :
437  Eps_ampl(eps_ampl), Pcos_initial(pcos_initial),
438  Pcos_duration(pcos_duration)
439 {
440  //-----------------------------------------------------------
441  // Create timesteppers
442  //-----------------------------------------------------------
443 
444  // Allocate the wall timestepper and add it to the problem's vector
445  // of timesteppers
446  Wall_time_stepper_pt = new Newmark<2>;
447  add_time_stepper_pt(Wall_time_stepper_pt);
448 
449  // Allocate the fluid timestepper and add it to the problem's Vector
450  // of timesteppers
451  Fluid_time_stepper_pt = new BDF<2>;
452  add_time_stepper_pt(Fluid_time_stepper_pt);
453 
454  //----------------------------------------------------------
455  // Create the wall mesh
456  //----------------------------------------------------------
457 
458  // Undeformed wall is an elliptical ring
459  Undef_geom_pt = new Ellipse(1.0,1.0);
460 
461  //Length of wall in Lagrangian coordinates
462  double L = 2.0*atan(1.0);
463 
464  //Now create the (Lagrangian!) mesh
465  Wall_mesh_pt = new
466  OneDLagrangianMesh<SOLID_ELEMENT>(N,L,Undef_geom_pt,Wall_time_stepper_pt);
467 
468  //----------------------------------------------------------
469  // Set the boundary conditions for wall mesh (problem)
470  //----------------------------------------------------------
471 
472  // Bottom boundary: (Boundary 0)
473  // No vertical displacement
474  Wall_mesh_pt->boundary_node_pt(0,0)->pin_position(1);
475  // Zero slope: Pin type 1 dof for displacement direction 0
476  Wall_mesh_pt->boundary_node_pt(0,0)->pin_position(1,0);
477 
478  // Top boundary: (Boundary 1)
479  // No horizontal displacement
480  Wall_mesh_pt->boundary_node_pt(1,0)->pin_position(0);
481  // Zero slope: Pin type 1 dof for displacement direction 1
482  Wall_mesh_pt->boundary_node_pt(1,0)->pin_position(1,1);
483 
484 
485  //-----------------------------------------------------------
486  // Create the fluid mesh:
487  //-----------------------------------------------------------
488 
489  // Fluid mesh is suspended from wall between the following Lagrangian
490  // coordinates:
491  double xi_lo=0.0;
492  double xi_hi=L;
493 
494  // Fractional position of dividing line for two outer blocks in mesh
495  double fract_mid=0.5;
496 
497  //Create a geometric object that represents the wall geometry from the
498  //wall mesh (one Lagrangian, two Eulerian coordinates).
499  MeshAsGeomObject *wall_mesh_as_geometric_object_pt
500  = new MeshAsGeomObject(Wall_mesh_pt);
501 
502  // Build fluid mesh using the wall mesh as a geometric object
503  Fluid_mesh_pt = new AlgebraicRefineableQuarterCircleSectorMesh<FLUID_ELEMENT >
504  (wall_mesh_as_geometric_object_pt,
505  xi_lo,fract_mid,xi_hi,Fluid_time_stepper_pt);
506 
507  // Set the error estimator
508  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
509  Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
510 
511  // Extract pointer to node at center of mesh
512  unsigned nnode=Fluid_mesh_pt->finite_element_pt(0)->nnode();
513  Veloc_trace_node_pt=Fluid_mesh_pt->finite_element_pt(0)->node_pt(nnode-1);
514 
515  //-------------------------------------------------------
516  // Set the fluid boundary conditions
517  //-------------------------------------------------------
518 
519  // Bottom boundary (boundary 0):
520  {
521  unsigned n_node = Fluid_mesh_pt->nboundary_node(0);
522  for (unsigned n=0;n<n_node;n++)
523  {
524  // Pin vertical velocity
525  Fluid_mesh_pt->boundary_node_pt(0,n)->pin(1);
526  }
527  }
528 
529  // Ring boundary (boundary 1):
530  // No slip; this also implies that the velocity needs
531  // to be updated in response to wall motion
532  {
533  unsigned n_node = Fluid_mesh_pt->nboundary_node(1);
534  for (unsigned n=0;n<n_node;n++)
535  {
536  // Which node are we dealing with?
537  Node* node_pt=Fluid_mesh_pt->boundary_node_pt(1,n);
538 
539  // Set auxiliary update function pointer
540  node_pt->set_auxiliary_node_update_fct_pt(
541  FSI_functions::apply_no_slip_on_moving_wall);
542 
543  // Pin both velocities
544  for(unsigned i=0;i<2;i++) {node_pt->pin(i);}
545  }
546  }
547 
548  // Left boundary (boundary 2):
549  {
550  unsigned n_node = Fluid_mesh_pt->nboundary_node(2);
551  for (unsigned n=0;n<n_node;n++)
552  {
553  // Pin horizontal velocity
554  Fluid_mesh_pt->boundary_node_pt(2,n)->pin(0);
555  }
556  }
557 
558 
559  //--------------------------------------------------------
560  // Add the submeshes and build global mesh
561  // -------------------------------------------------------
562 
563  // Wall mesh
564  add_sub_mesh(Wall_mesh_pt);
565 
566  //Fluid mesh
567  add_sub_mesh(Fluid_mesh_pt);
568 
569  // Combine all submeshes into a single Mesh
570  build_global_mesh();
571 
572 
573  //----------------------------------------------------------
574  // Finish problem setup
575  // ---------------------------------------------------------
576 
577  //Find number of elements in fluid mesh
578  unsigned n_element = Fluid_mesh_pt->nelement();
579 
580  // Loop over the fluid elements to set up element-specific
581  // things that cannot be handled by constructor
582  for(unsigned e=0;e<n_element;e++)
583  {
584  // Upcast from FiniteElement to the present element
585  FLUID_ELEMENT *el_pt
586  = dynamic_cast<FLUID_ELEMENT*>(Fluid_mesh_pt->element_pt(e));
587 
588  //Set the Reynolds number, etc
589  el_pt->re_pt() = &Global_Physical_Variables::Re;
590  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
591 
592 
593  el_pt->evaluate_shape_derivs_by_direct_fd();
594 
595 // el_pt->evaluate_shape_derivs_by_chain_rule();
596 // el_pt->enable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
597 
598 // if (e==0)
599 // {
600 // el_pt->disable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
601 // }
602 // else
603 // {
604 // el_pt->enable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
605 // }
606 
607  //el_pt->evaluate_shape_derivs_by_direct_fd();
608 
609  }
610 
611 
612  //Loop over the solid elements to set physical parameters etc.
613  unsigned n_wall_element = Wall_mesh_pt->nelement();
614  for(unsigned e=0;e<n_wall_element;e++)
615  {
616  //Cast to proper element type
617  SOLID_ELEMENT *el_pt = dynamic_cast<SOLID_ELEMENT*>(
618  Wall_mesh_pt->element_pt(e));
619 
620  //Set physical parameters for each element:
621  el_pt->h_pt() = &Global_Physical_Variables::H;
622  el_pt->lambda_sq_pt() = &Global_Physical_Variables::Lambda_sq;
623 
624  //Function that specifies the external load Vector
625  el_pt->load_vector_fct_pt() = &Global_Physical_Variables::pcos_load;
626 
627  // Function that specifies the load ratios
628  el_pt->q_pt() = &Global_Physical_Variables::Q;
629 
630  //Assign the undeformed beam shape
631  el_pt->undeformed_beam_pt() = Undef_geom_pt;
632  }
633 
634  // Establish control displacment: (even if no displacement control is applied
635  // we still want to doc the displacement at the same point)
636 
637  // Choose element: (This is the last one)
638  Doc_displacement_elem_pt = dynamic_cast<SOLID_ELEMENT*>(
639  Wall_mesh_pt->element_pt(n_wall_element-1));
640 
641  // Setup fsi: Work out which fluid dofs affect the wall elements
642  // the correspondance between wall dofs and fluid elements is handled
643  // during the remeshing, but the "reverse" association must be done
644  // separately.
645  // We pass the boundary between the fluid and solid meshes and pointers
646  // to the meshes. The interaction boundary is boundary 1 of the
647  // 2D fluid mesh.
648  FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
649  (this,1,Fluid_mesh_pt,Wall_mesh_pt);
650 
651  // Do equation numbering
652  cout << "# of dofs " << assign_eqn_numbers() << std::endl;
653 
654 }
655 
656 
657 //=========================================================================
658 /// Solver loop to perform unsteady run
659 //=========================================================================
661 {
662  // Setup documentation
663  //---------------------------------------------------------------
664 
665  /// Label for output
666  DocInfo doc_info;
667 
668  // Output directory
669  doc_info.set_directory("RESLT");
670 
671  // Step number
672  doc_info.number()=0;
673 
674  //Open a trace file
675  ofstream trace_file("RESLT/trace_ring.dat");
676 
677  // Write header for trace file
678  trace_file << "VARIABLES=\"time\",\"V_c_t_r_l\"";
679  trace_file << ",\"x<sub>1</sub><sup>(track)</sup>\"";
680  trace_file << ",\"x<sub>2</sub><sup>(track)</sup>\"";
681  trace_file << ",\"u<sub>1</sub><sup>(track)</sup>\"";
682  trace_file << ",\"u<sub>2</sub><sup>(track)</sup>\"";
683  trace_file << ",\"N<sub>element</sub>\"";
684  trace_file << ",\"N<sub>dof</sub>\"";
685  trace_file << ",\"# of under-refined elements\"";
686  trace_file << ",\"max. error\"";
687  trace_file << ",\"min. error\"";
688  trace_file << ",\"max. permitted error\"";
689  trace_file << ",\"min. permitted error\"";
690  trace_file << ",\"max. permitted # of unrefined elements\"";
691  trace_file << ",\"doc number\"";
692  trace_file << std::endl;
693 
694 
695  // Initialise timestepping
696  // -------------------------------------------------------------
697 
698  // Number of steps
699  unsigned nstep=300;
700 
701  // Nontrivial command line input: validation: only do three steps
702  if (CommandLineArgs::Argc>1)
703  {
704  nstep=1;
705  cout << "Only doing nstep steps for validation: " << nstep << std::endl;
706  }
707 
708  // Set initial timestep
709  double dt=0.004;
710 
711  // Set initial value for dt -- also assigns weights to the timesteppers
712  initialise_dt(dt);
713 
714  // Set physical parameters
715  //---------------------------------------------------------
716  using namespace Global_Physical_Variables;
717 
718  // Set Womersley number
719  Alpha_sq=100.0; // 50.0; // ADJUST
720 
721  // Set density ratio
722  Density_ratio=10.0; // 0.0; ADJUST
723 
724  // Wall thickness
725  H=1.0/20.0;
726 
727  // Set external pressure
728  Pext=0.0;
729 
730  /// Perturbation pressure
732 
733  // Assign/doc corresponding computational parameters
734  set_params();
735 
736 
737  // Refine uniformly and assign initial conditions
738  //--------------------------------------------------------------
739 
740  // Refine the problem uniformly
741  refine_uniformly();
742  refine_uniformly();
743 
744  // This sets up the solution at the initial time
746 
747  // Set targets for spatial adptivity
748  //---------------------------------------------------------------
749 
750  // Max. and min. error for adaptive refinement/unrefinement
751  Fluid_mesh_pt->max_permitted_error()=1.0e-2;
752  Fluid_mesh_pt->min_permitted_error()=1.0e-3;
753 
754  // Don't allow refinement to drop under given level
755  Fluid_mesh_pt->min_refinement_level()=2;
756 
757  // Don't allow refinement beyond given level
758  Fluid_mesh_pt->max_refinement_level()=6;
759 
760  // Don't bother adapting the mesh if no refinement is required
761  // and if less than ... elements are to be merged.
762  Fluid_mesh_pt->max_keep_unrefined()=20;
763 
764  // Doc refinement targets
765  Fluid_mesh_pt->doc_adaptivity_targets(cout);
766 
767 
768  // Do the timestepping
769  //----------------------------------------------------------------
770 
771  // Reset initial conditions after refinement for first step only
772  bool first=true;
773 
774  //Output initial data
775  doc_solution(0,doc_info,trace_file);
776 
777 
778 // {
779 // unsigned nel=Fluid_mesh_pt->nelement();
780 // for (unsigned e=0;e<nel;e++)
781 // {
782 // std::cout << "\n\nEl: " << e << std::endl << std::endl;
783 // FiniteElement* el_pt=Fluid_mesh_pt->finite_element_pt(e);
784 // unsigned n_dof=el_pt->ndof();
785 // Vector<double> residuals(n_dof);
786 // DenseDoubleMatrix jacobian(n_dof,n_dof);
787 // el_pt->get_jacobian(residuals,jacobian);
788 // }
789 // exit(0);
790 // }
791 
792  // Time integration loop
793  for(unsigned i=1;i<=nstep;i++)
794  {
795  // Switch doc off during solve
796  doc_info.disable_doc();
797 
798  // Solve
799  unsigned max_adapt=1;
800  unsteady_newton_solve(dt,max_adapt,first);
801 
802  // Now we've done the first step
803  first=false;
804 
805  // Doc solution
806  doc_solution(i,doc_info,trace_file);
807 
808  /// Switch off perturbation pressure
809  if (time_pt()->time()>Pcos_duration) {Pcos=0.0;}
810 
811  }
812 
813 }
814 
815 
816 //=====================================================================
817 /// Driver for fsi ring test problem
818 //=====================================================================
819 int main(int argc, char* argv[])
820 {
821 
822  // Store command line arguments
823  CommandLineArgs::setup(argc,argv);
824 
825  // Number of elements
826  unsigned nelem = 13;
827 
828  /// Perturbation pressure
829  double pcos_initial=1.0e-6; // ADJUST
830 
831  /// Duration of initial pcos perturbation
832  double pcos_duration=0.3; // ADJUST
833 
834  /// Amplitude of initial deformation
835  double eps_ampl=0.0; // ADJUST
836 
837  //Set up the problem
838  FSIRingProblem problem(nelem,eps_ampl,pcos_initial,pcos_duration);
839 
840  // Do parameter study
841  problem.dynamic_run();
842 
843 }
844 
845 
846 
847 
848 
849 
void doc_solution(const unsigned &i, DocInfo &doc_info, ofstream &trace_file)
Doc solution: Pass number of timestep, i (we append to tracefile after every timestep but do a full d...
double Density_ratio
Density ratio of the solid and the fluid.
Definition: fsi_osc_ring.cc:59
void set_params()
Set the parameters that are used in the code as a function of the Womersley number, the density ratio and H.
Definition: fsi_osc_ring.cc:91
double Lambda_sq
Timescale ratio (non-dimensation density)
Definition: fsi_osc_ring.cc:84
double Pcos
Perturbation pressure.
Definition: fsi_osc_ring.cc:71
AlgebraicRefineableQuarterCircleSectorMesh< FLUID_ELEMENT > * Fluid_mesh_pt
Pointer to fluid mesh.
void set_fluid_initial_condition()
Setup initial condition for fluid.
double Pcos_initial
Initial pcos.
double Q
Stress ratio.
Definition: fsi_osc_ring.cc:87
void dynamic_run()
Do dynamic run.
double ReSt
Reynolds x Strouhal number.
Definition: fsi_osc_ring.cc:81
void set_initial_condition()
Setup initial condition for both domains.
double Eps_ampl
Amplitude of initial deformation.
double Pcos_duration
Duration of initial pcos.
double Pext
External Pressure.
Definition: fsi_osc_ring.cc:62
double Re
Reynolds number.
Definition: fsi_osc_ring.cc:78
BDF< 2 > * Fluid_time_stepper_pt
Pointer to fluid timestepper.
double H
Nondimensional thickness of the beam.
Definition: fsi_osc_ring.cc:68
GeomObject * Undef_geom_pt
Pointer to geometric object that represents the undeformed wall shape.
SOLID_ELEMENT * Doc_displacement_elem_pt
Element used for documenting displacement.
Newmark< 2 > * Wall_time_stepper_pt
Pointer to wall timestepper.
FSIHermiteBeamElement SOLID_ELEMENT
Typedef to specify the solid element used.
void actions_after_adapt()
Update the problem specs after adaptation:
double Alpha_sq
Square of Womersly number (a frequency parameter)
Definition: fsi_osc_ring.cc:56
void actions_after_newton_solve()
Update after solve (empty)
int main(int argc, char *argv[])
Driver for fsi ring test problem.
FSIRingProblem(const unsigned &nelement_wall, const double &eps_ampl, const double &pcos_initial, const double &pcos_duration)
void pcos_load(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
void set_wall_initial_condition()
Setup initial condition for wall.
AlgebraicElement< RefineableQCrouzeixRaviartElement< 2 > > FLUID_ELEMENT
There are very few element types that will work for this problem. Rather than passing the element typ...
OneDLagrangianMesh< SOLID_ELEMENT > * Wall_mesh_pt
Pointer to wall mesh.
double Nu
Poisson ratio.
Definition: fsi_osc_ring.cc:65
void actions_before_newton_solve()
Update before solve (empty)
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.