unsteady_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 code for a large-amplitude ring oscillation
31 
32 //Generic stuff
33 #include "generic.h"
34 
35 // The beam equations
36 #include "beam.h"
37 
38 // The mesh
39 #include "meshes/one_d_lagrangian_mesh.h"
40 
41 using namespace std;
42 
43 using namespace oomph;
44 
45 
46 ////////////////////////////////////////////////////////////////////////
47 ////////////////////////////////////////////////////////////////////////
48 ////////////////////////////////////////////////////////////////////////
49 
50 
51 //====start_of_namespace============================
52 /// Namespace for global parameters
53 //==================================================
54 namespace Global_Physical_Variables
55 {
56 
57  /// Perturbation pressure
58  double Pcos;
59 
60  /// Duration of transient load
61  double T_kick;
62 
63  /// Load function: Perturbation pressure to force non-axisymmetric deformation
64  void press_load(const Vector<double>& xi,
65  const Vector<double> &x,
66  const Vector<double>& N,
67  Vector<double>& load)
68  {
69  for(unsigned i=0;i<2;i++)
70  {
71  load[i] = -Pcos*cos(2.0*xi[0])*N[i];
72  }
73  } //end of load
74 
75 
76  /// Scaling factor for wall thickness (to be used in an exercise)
77  double Alpha=1.0;
78 
79  /// Wall thickness -- 1/20 for default value of scaling factor
80  double H=Alpha*1.0/20.0;
81 
82  /// \short Square of timescale ratio (i.e. non-dimensional density)
83  /// -- 1.0 for default value of scaling factor
84  double Lambda_sq=pow(Alpha,2);
85 
86 } // end of namespace
87 
88 
89 
90 
91 
92 /////////////////////////////////////////////////////////////////////
93 /////////////////////////////////////////////////////////////////////
94 /////////////////////////////////////////////////////////////////////
95 
96 
97 
98 //===start_of_problem_class=============================================
99 /// Ring problem
100 //======================================================================
101 template<class ELEMENT, class TIMESTEPPER>
102 class ElasticRingProblem : public Problem
103 {
104 
105 public:
106 
107  /// Constructor: Number of elements
108  ElasticRingProblem(const unsigned& n_element);
109 
110  /// Access function for the specific mesh
111  OneDLagrangianMesh<ELEMENT>* mesh_pt()
112  {
113  return dynamic_cast<OneDLagrangianMesh<ELEMENT>*>(Problem::mesh_pt());
114  }
115 
116  /// Update function is empty
118 
119  /// Update function is empty
121 
122  /// \short Setup initial conditions
123  void set_initial_conditions();
124 
125  /// Doc solution
126  void doc_solution(DocInfo& doc_info);
127 
128  /// Do unsteady run
129  void unsteady_run();
130 
131  /// \short Dump problem-specific parameter values, then dump
132  /// generic problem data.
133  void dump_it(ofstream& dump_file);
134 
135  /// \short Read problem-specific parameter values, then recover
136  /// generic problem data.
137  void restart(ifstream& restart_file);
138 
139 
140 private:
141 
142  /// Trace file for recording control data
143  ofstream Trace_file;
144 
145  /// Flag for validation run: Default: 0 = no validation run
147 
148  /// Restart flag specified via command line?
150 
151 }; // end of problem class
152 
153 
154 
155 
156 
157 //===start_of_constructor===============================================
158 /// Constructor for elastic ring problem
159 //======================================================================
160 template<class ELEMENT,class TIMESTEPPER>
162 (const unsigned& n_element) :
163  Validation_run_flag(0), //default: false
164  Restart_flag(false)
165 {
166 
167  // Create the timestepper and add it to the Problem's collection of
168  // timesteppers -- this creates the Problem's Time object.
169  add_time_stepper_pt(new TIMESTEPPER());
170 
171  // Undeformed beam is an ellipse with unit axes
172  GeomObject* undef_geom_pt=new Ellipse(1.0,1.0);
173 
174  //Length of domain
175  double length = MathematicalConstants::Pi/2.0;
176 
177  //Now create the (Lagrangian!) mesh
178  Problem::mesh_pt() = new OneDLagrangianMesh<ELEMENT>(
179  n_element,length,undef_geom_pt,Problem::time_stepper_pt());
180 
181  // Boundary condition:
182 
183  // Bottom:
184  unsigned ibound=0;
185  // No vertical displacement
186  mesh_pt()->boundary_node_pt(ibound,0)->pin_position(1);
187  // Zero slope: Pin type 1 (slope) dof for displacement direction 0
188  mesh_pt()->boundary_node_pt(ibound,0)->pin_position(1,0);
189 
190  // Top:
191  ibound=1;
192  // No horizontal displacement
193  mesh_pt()->boundary_node_pt(ibound,0)->pin_position(0);
194  // Zero slope: Pin type 1 (slope) dof for displacement direction 1
195  mesh_pt()->boundary_node_pt(ibound,0)->pin_position(1,1);
196 
197  // Complete build of all elements so they are fully functional
198  // -----------------------------------------------------------
199 
200  //Loop over the elements to set physical parameters etc.
201  for(unsigned i=0;i<n_element;i++)
202  {
203  // Cast to proper element type
204  ELEMENT *elem_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
205 
206  // Pass pointer to square of timescale ratio (non-dimensional density)
207  elem_pt->lambda_sq_pt() = &Global_Physical_Variables::Lambda_sq;
208 
209  // Pass pointer to non-dimensional wall thickness
210  elem_pt->h_pt() = &Global_Physical_Variables::H;
211 
212  // Function that specifies load vector
213  elem_pt->load_vector_fct_pt() = &Global_Physical_Variables::press_load;
214 
215  // Assign the undeformed surface
216  elem_pt->undeformed_beam_pt() = undef_geom_pt;
217  }
218 
219  // Do equation numbering
220  cout << "# of dofs " << assign_eqn_numbers() << std::endl;
221 
222 } // end of constructor
223 
224 
225 
226 
227 //====start_of_doc_solution===============================================
228 /// Document solution
229 //========================================================================
230 template<class ELEMENT, class TIMESTEPPER>
232  DocInfo& doc_info)
233 {
234 
235  cout << "Doc-ing step " << doc_info.number()
236  << " for time " << time_stepper_pt()->time_pt()->time() << std::endl;
237 
238 
239  // Loop over all elements to get global kinetic and potential energy
240  unsigned n_elem=mesh_pt()->nelement();
241  double global_kin=0;
242  double global_pot=0;
243  double pot,kin;
244  for (unsigned ielem=0;ielem<n_elem;ielem++)
245  {
246  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(ielem))->get_energy(pot,kin);
247  global_kin+=kin;
248  global_pot+=pot;
249  }
250 
251 
252  // Get pointer to last element to document displacement
253  FiniteElement* trace_elem_pt=mesh_pt()->finite_element_pt(n_elem-1);
254 
255  // Vector of local coordinates at control point
256  Vector<double> s_trace(1);
257  s_trace[0]=1.0;
258 
259  // Write trace file: Time, control position, energies
260  Trace_file << time_pt()->time() << " "
261  << trace_elem_pt->interpolated_x(s_trace,1)
262  << " " << global_pot << " " << global_kin
263  << " " << global_pot + global_kin
264  << std::endl; // end of output to trace file
265 
266  ofstream some_file;
267  char filename[100];
268 
269  // Number of plot points
270  unsigned npts=5;
271 
272  // Output solution
273  sprintf(filename,"%s/ring%i.dat",doc_info.directory().c_str(),
274  doc_info.number());
275  some_file.open(filename);
276  mesh_pt()->output(some_file,npts);
277 
278  // Write file as a tecplot text object
279  some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
280  << time_pt()->time() << "\"";
281 
282  // ...and draw a horizontal line whose length is proportional
283  // to the elapsed time
284  some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
285  some_file << "1" << std::endl;
286  some_file << "2" << std::endl;
287  some_file << " 0 0" << std::endl;
288  some_file << time_pt()->time()*20.0 << " 0" << std::endl;
289  some_file.close();
290 
291 
292  // Loop over all elements do dump out previous solutions
293  unsigned nsteps=time_stepper_pt()->nprev_values();
294  for (unsigned t=0;t<=nsteps;t++)
295  {
296  sprintf(filename,"%s/ring%i-%i.dat",doc_info.directory().c_str(),
297  doc_info.number(),t);
298  some_file.open(filename);
299  unsigned n_elem=mesh_pt()->nelement();
300  for (unsigned ielem=0;ielem<n_elem;ielem++)
301  {
302  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(ielem))->
303  output(t,some_file,npts);
304  }
305  some_file.close();
306  } // end of output of previous solutions
307 
308 
309  // Write restart file
310  sprintf(filename,"%s/ring_restart%i.dat",doc_info.directory().c_str(),
311  doc_info.number());
312  some_file.open(filename);
313  dump_it(some_file);
314  some_file.close();
315 
316 
317 } // end of doc solution
318 
319 
320 
321 
322 
323 
324 //===start_of_dump_it====================================================
325 /// Dump problem-specific parameter values, then dump
326 /// generic problem data.
327 //=======================================================================
328 template<class ELEMENT,class TIMESTEPPER>
330 
331 {
332 
333  // Write Pcos
334  dump_file << Global_Physical_Variables::Pcos << " # Pcos" << std::endl;
335 
336  // Write validation run flag
337  dump_file << Validation_run_flag
338  << " # Validation run flag" << std::endl;
339 
340  // Call generic dump()
341  Problem::dump(dump_file);
342 
343 } // end of dump it
344 
345 
346 //==start_of_restart=====================================================
347 /// Read problem-specific parameter values, then recover
348 /// generic problem data.
349 //=======================================================================
350 template<class ELEMENT,class TIMESTEPPER>
352 {
353 
354  string input_string;
355 
356  // Read line up to termination sign
357  getline(restart_file,input_string,'#');
358  // Ignore rest of line
359  restart_file.ignore(80,'\n');
360  // Read in pcos
361  Global_Physical_Variables::Pcos=atof(input_string.c_str());
362 
363  // Read line up to termination sign
364  getline(restart_file,input_string,'#');
365  // Ignore rest of line
366  restart_file.ignore(80,'\n');
367  // Read in Long run flag
368  Validation_run_flag=
369  unsigned(atof(input_string.c_str()));
370 
371  // Call generic read()
372  Problem::read(restart_file);
373 
374 } // end of restart
375 
376 
377 
378 
379 //=====start_of_set_ic=====================================================
380 /// Setup initial conditions -- either restart from solution
381 /// specified via command line or impulsive start.
382 //=========================================================================
383 template<class ELEMENT,class TIMESTEPPER>
385 {
386 
387  // No restart file --> impulsive start from initial configuration
388  // assigned in the Lagrangian mesh.
389  if (!Restart_flag)
390  {
391  // Set initial timestep
392  double dt=1.0;
393 
394  // Assign impulsive start
395  assign_initial_values_impulsive(dt);
396  }
397  // Restart file specified via command line
398  else
399  {
400  // Try to open restart file
401  ifstream* restart_file_pt=
402  new ifstream(CommandLineArgs::Argv[2],ios_base::in);
403  if (restart_file_pt!=0)
404  {
405  cout << "Have opened " << CommandLineArgs::Argv[2] <<
406  " for restart. " << std::endl;
407  }
408  else
409  {
410  std::ostringstream error_stream;
411  error_stream << "ERROR while trying to open "
412  << CommandLineArgs::Argv[2]
413  << " for restart." << std::endl;
414 
415  throw OomphLibError(error_stream.str(),
416  OOMPH_CURRENT_FUNCTION,
417  OOMPH_EXCEPTION_LOCATION);
418  }
419 
420  // Read restart data:
421  restart(*restart_file_pt);
422 
423  }
424 
425 } // end of set ic
426 
427 
428 
429 //===start_of_unsteady_run=================================================
430 /// Solver loop to perform unsteady run.
431 //=========================================================================
432 template<class ELEMENT,class TIMESTEPPER>
434 {
435 
436  // Convert command line arguments (if any) into flags:
437  //----------------------------------------------------
438 
439  if (CommandLineArgs::Argc<2)
440  {
441  cout << "No command line arguments -- using defaults."
442  << std::endl;
443  }
444  else if (CommandLineArgs::Argc==2)
445  {
446  // Flag for validation run
447  Validation_run_flag=atoi(CommandLineArgs::Argv[1]);
448  }
449  else if (CommandLineArgs::Argc==3)
450  {
451  // Flag for validation run
452  Validation_run_flag=atoi(CommandLineArgs::Argv[1]);
453 
454  // Second argument is restart file. If it's there we're performing
455  // a restart
456  Restart_flag=true;
457  }
458  else
459  {
460  std::string error_message =
461  "Wrong number of command line arguments. Specify two or fewer.\n";
462  error_message += "Arg1: Validation_run_flag [0/1] for [false/true]\n";
463  error_message += "Arg2: Name of restart_file (optional)\n";
464  error_message += "No arguments specified --> default run\n";
465 
466  throw OomphLibError(error_message,
467  OOMPH_CURRENT_FUNCTION,
468  OOMPH_EXCEPTION_LOCATION);
469  }
470 
471 
472  // Label for output
473  DocInfo doc_info;
474 
475  // Output directory
476  doc_info.set_directory("RESLT");
477 
478  // Step number
479  doc_info.number()=0;
480 
481  // Set up trace file
482  char filename[100];
483  sprintf(filename,"%s/trace_ring.dat",doc_info.directory().c_str());
484  Trace_file.open(filename);
485 
486  Trace_file << "VARIABLES=\"time\",\"R<sub>ctrl</sub>\",\"E<sub>pot</sub>\"";
487  Trace_file << ",\"E<sub>kin</sub>\",\"E<sub>kin</sub>+E<sub>pot</sub>\""
488  << std::endl;
489 
490  // Perturbation pressure -- incl. scaling factor (Alpha=1.0 by default)
493 
494  // Duration of transient load
496 
497  // Number of steps
498  unsigned nstep=600;
499  if (Validation_run_flag==1) {nstep=10;}
500 
501  // Setup initial condition (either restart or impulsive start)
502  set_initial_conditions();
503 
504  // Extract initial timestep as set up in set_initial_conditions()
505  double dt=time_pt()->dt();
506 
507  // Output initial data
508  doc_solution(doc_info);
509 
510  // If the run is restarted, dt() contains the size of the timestep
511  // that was taken to reach the dumped solution. In a non-restarted run
512  // the next step would have been taken with a slightly smaller
513  // timestep (see below). For a restarted run we therefore adjust
514  // the next timestep accordingly.
515  if (Restart_flag&&(time_pt()->time()>10.0)&&(time_pt()->time()<100.0))
516  {
517  dt=0.995*dt;
518  }
519 
520  // Time integration loop
521  for(unsigned i=1;i<=nstep;i++)
522  {
523  // Switch off perturbation pressure
524  if (time_pt()->time()>Global_Physical_Variables::T_kick)
525  {
526  /// Perturbation pressure
528  }
529 
530  // Solve
531  unsteady_newton_solve(dt);
532 
533  // Doc solution
534  doc_info.number()++;
535  doc_solution(doc_info);
536 
537  // Reduce timestep for the next solve
538  if ((time_pt()->time()>10.0)&&(time_pt()->time()<100.0))
539  {
540  dt=0.995*dt;
541  }
542 
543  }
544 
545 } // end of unsteady run
546 
547 
548 
549 
550 //===start_of_main=====================================================
551 /// Driver for oscillating ring problem
552 //=====================================================================
553 int main(int argc, char* argv[])
554 {
555 
556  // Store command line arguments
557  CommandLineArgs::setup(argc,argv);
558 
559  // Number of elements
560  unsigned nelem = 13;
561 
562  //Set up the problem
564 
565  // Do unsteady run
566  problem.unsteady_run();
567 
568 } // end of main
569 
570 
571 
572 
573 
574 
575 
576 
double Lambda_sq
Square of timescale ratio (i.e. non-dimensional density) – 1.0 for default value of scaling factor...
double Pcos
Perturbation pressure.
void set_initial_conditions()
Setup initial conditions.
void unsteady_run()
Do unsteady run.
unsigned Validation_run_flag
Flag for validation run: Default: 0 = no validation run.
double T_kick
Duration of transient load.
void press_load(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Load function: Perturbation pressure to force non-axisymmetric deformation.
int main(int argc, char *argv[])
Driver for oscillating ring problem.
void doc_solution(DocInfo &doc_info)
Doc solution.
OneDLagrangianMesh< ELEMENT > * mesh_pt()
Access function for the specific mesh.
void actions_before_newton_solve()
Update function is empty.
void dump_it(ofstream &dump_file)
Dump problem-specific parameter values, then dump generic problem data.
bool Restart_flag
Restart flag specified via command line?
double H
Wall thickness – 1/20 for default value of scaling factor.
void restart(ifstream &restart_file)
Read problem-specific parameter values, then recover generic problem data.
void actions_after_newton_solve()
Update function is empty.
double Alpha
Scaling factor for wall thickness (to be used in an exercise)
ElasticRingProblem(const unsigned &N, const double &L)
Constructor: Number of elements, length of domain, flag for setting Newmark IC directly or consistent...