bretherton.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 // The 2D Bretherton problem
31 #include<algorithm>
32 
33 // The oomphlib headers
34 #include "generic.h"
35 #include "navier_stokes.h"
36 #include "fluid_interface.h"
37 
38 // The mesh
39 #include "meshes/bretherton_spine_mesh.h"
40 
41 using namespace std;
42 
43 using namespace oomph;
44 
45 //======================================================================
46 /// Namepspace for global parameters
47 //======================================================================
48 namespace Global_Physical_Variables
49 {
50 
51  /// Reynolds number
52  double Re;
53 
54  /// Womersley = Reynolds times Strouhal
55  double ReSt;
56 
57  /// Product of Reynolds and Froude number
58  double ReInvFr;
59 
60  /// Capillary number
61  double Ca;
62 
63  /// Direction of gravity
64  Vector<double> G(2);
65 
66  /// Pointer to film thickness at outflow on the lower wall
67  double* H_lo_pt;
68 
69  /// Pointer to film thickness at outflow on the upper wall
70  double* H_up_pt;
71 
72  /// Pointer to y-position at inflow on the lower wall
73  double* Y_lo_pt;
74 
75  /// Pointer to y-position at inflow on the upper wall
76  double* Y_up_pt;
77 
78  /// Set inflow velocity, based on spine heights at outflow
79  /// and channel width at inflow
80  void inflow(const Vector<double>& x, Vector<double>& veloc)
81  {
82 #ifdef PARANOID
83  std::ostringstream error_stream;
84  bool throw_error=false;
85  if (H_lo_pt==0)
86  {
87  error_stream << "You must set H_lo_pt\n";
88  throw_error = true;
89  }
90  if (H_up_pt==0)
91  {
92  error_stream << "You must set H_up_pt\n";
93  throw_error = true;
94  }
95  if (Y_lo_pt==0)
96  {
97  error_stream << "You must set Y_lo_pt\n";
98  throw_error = true;
99  }
100  if (Y_up_pt==0)
101  {
102  error_stream << "You must set Y_up_pt\n";
103  throw_error = true;
104  }
105 
106  //If one of the errors has occured
107  if(throw_error)
108  {
109  throw OomphLibError(error_stream.str(),
110  OOMPH_CURRENT_FUNCTION,
111  OOMPH_EXCEPTION_LOCATION);
112  }
113 #endif
114 
115 
116  // Get average film thickness far ahead
117  double h_av=0.5*(*H_lo_pt+*H_up_pt);
118 
119  // Get position of upper and lower wall at inflow
120  double y_up=*Y_up_pt;
121  double y_lo=*Y_lo_pt;
122 
123  // Constant in velocity profile
124  double C =6.0*(2.0*h_av+y_lo-y_up)/
125  (y_up*y_up*y_up-y_lo*y_lo*y_lo-h_av*y_up*
126  y_up*y_up+h_av*y_lo*y_lo*y_lo-3.0*y_lo*y_up*y_up+
127  3.0*y_lo*y_lo*y_up+3.0*y_lo*y_up*y_up*h_av-3.0*y_lo*y_lo*y_up*h_av);
128 
129  // y coordinate
130  double y=x[1];
131 
132  // Parallel, parabolic inflow
133  veloc[0]=-1.0+C*(1.0-h_av)*((y_lo-y)*(y_up-y));
134  veloc[1]=0.0;
135  }
136 
137 }
138 
139 
140 
141 
142 //======================================================================
143 /// "Bretherton element" is a fluid element (of type ELEMENT) for which
144 /// the (inflow) velocity at those nodes that are located on a
145 /// specified Mesh boundary is prescribed by Dirichlet boundary
146 /// conditions. The key is that prescribed velocity profile can be
147 /// a function of some external Data -- this dependency must
148 /// be taken into account when computing the element's Jacobian
149 /// matrix.
150 ///
151 /// This element type is useful, for instance, in the Bretherton
152 /// problem, where the parabolic "inflow" profile is a function
153 /// of the film thickness (represented by Spine heights) at the
154 /// "outflow".
155 //======================================================================
156 template<class ELEMENT>
157 class BrethertonElement : public ELEMENT
158 {
159 
160 public:
161 
162  /// \short Typedef for pointer (global) function that specifies the
163  /// the inflow
164  typedef void (*InflowFctPt)(const Vector<double>& x, Vector<double>& veloc);
165 
166  /// Constructor: Call the constructor of the underlying element
167  BrethertonElement() : ELEMENT() {}
168 
169 
170  /// \short Activate the dependency of the "inflow" on the external
171  /// data. Pass the vector of pointers to the external Data that affects
172  /// the inflow, the id of the boundary on which the inflow
173  /// condition is to be applied and the function pointer to
174  /// to the global function that defines the inflow
176  const Vector<Data*>& inflow_ext_data,
177  const unsigned& inflow_boundary,
178  InflowFctPt inflow_fct_pt)
179  {
180  // Copy data that affects the inflow
181  unsigned n_ext=inflow_ext_data.size();
182  Inflow_ext_data.resize(n_ext);
183  for (unsigned i=0;i<n_ext;i++)
184  {
185  Inflow_ext_data[i]=inflow_ext_data[i];
186  }
187  // Set inflow boundary
188  Inflow_boundary=inflow_boundary;
189  // Set fct pointer to inflow condition
190  Inflow_fct_pt=inflow_fct_pt;
191  }
192 
193 
194  /// short Overload assign local equation numbers: Add the dependency
195  /// on the external Data that affects the inflow profile
196  void assign_local_eqn_numbers(const bool &store_local_dof_pt)
197  {
198  // Call the element's local equation numbering procedure first
199  ELEMENT::assign_local_eqn_numbers(store_local_dof_pt);
200 
201  // Now add the equation numbers for the Data that affects the inflow
202  // profile
203 
204  // Number of local equations so far
205  unsigned local_eqn_count = this->ndof();
206 
207  // Find out max. number of values stored at all Data values
208  // that affect the inflow
209  unsigned max_nvalue=0;
210  unsigned n_ext=Inflow_ext_data.size();
211  for (unsigned i=0;i<n_ext;i++)
212  {
213  // The external Data:
214  Data* data_pt=Inflow_ext_data[i];
215  // Number of values
216  unsigned n_val=data_pt->nvalue();
217  if (n_val>max_nvalue) max_nvalue=n_val;
218  }
219 
220  // Allocate sufficient storage
221  Inflow_ext_data_eqn.resize(n_ext,max_nvalue);
222 
223  //A local queue to store the global equation numbers
224  std::deque<unsigned long> global_eqn_number_queue;
225 
226  // Loop over external Data that affect the "inflow"
227  for (unsigned i=0;i<n_ext;i++)
228  {
229  // The external Data:
230  Data* data_pt=Inflow_ext_data[i];
231 
232  // Loop over number of values:
233  unsigned n_val=data_pt->nvalue();
234  for (unsigned ival=0;ival<n_val;ival++)
235  {
236 
237  // Is it free or pinned?
238  long eqn_number=data_pt->eqn_number(ival);
239  if (eqn_number>=0)
240  {
241  // Add global equation number to local storage
242  global_eqn_number_queue.push_back(eqn_number);
243  // Store local equation number associated with this external dof
244  Inflow_ext_data_eqn(i,ival)=local_eqn_count;
245  // Increment counter for local dofs
246  local_eqn_count++;
247  }
248  else
249  {
250  Inflow_ext_data_eqn(i,ival)=-1;
251  }
252  }
253  }
254 
255  //Now add our global equations numbers to the internal element storage
256  this->add_global_eqn_numbers(global_eqn_number_queue,
257  GeneralisedElement::Dof_pt_deque);
258  }
259 
260 
261  /// Overloaded Jacobian computation: Computes the Jacobian
262  /// of the underlying element and then adds the FD
263  /// operations to evaluate the derivatives w.r.t. the Data values
264  /// that affect the inflow.
265  void get_jacobian(Vector<double>& residuals,
266  DenseMatrix<double>& jacobian)
267  {
268  // Loop over Data values that affect the inflow
269  unsigned n_ext=Inflow_ext_data.size();
270 
271  // Call "normal" jacobian first.
272  ELEMENT::get_jacobian(residuals,jacobian);
273 
274  if (n_ext==0) return;
275 
276  // Get ready for FD operations
277  Vector<double> residuals_plus(residuals.size());
278  double fd_step=1.0e-8;
279 
280  // Loop over Data values that affect the inflow
281  for (unsigned i=0;i<n_ext;i++)
282  {
283  // Get Data item
284  Data* data_pt=Inflow_ext_data[i];
285 
286  // Loop over values
287  unsigned n_val=data_pt->nvalue();
288  for (unsigned ival=0;ival<n_val;ival++)
289  {
290  // Local equation number
291  int local_eqn=Inflow_ext_data_eqn(i,ival);
292 
293  // Dof or pinned?
294  if (local_eqn>=0)
295  {
296  //get pointer to the data value
297  double *value_pt = data_pt->value_pt(ival);
298 
299  //Backup Data value
300  double backup = *value_pt;
301 
302  // Do FD step
303  *value_pt += fd_step;
304 
305  // Re-assign the inflow velocities for nodes in this element
306  reassign_inflow();
307 
308 
309  // Fill in the relevant column in the Jacobian matrix
310  unsigned n_dof = this->ndof();
311  //Zero the residuals
312  for(unsigned idof=0;idof<n_dof;idof++) {residuals_plus[idof] = 0.0;}
313  // Re-compute the element residual
314  this->get_residuals(residuals_plus);
315 
316  for(unsigned idof=0;idof<n_dof;idof++)
317  {
318  jacobian(idof,local_eqn)=
319  (residuals_plus[idof]-residuals[idof])/fd_step;
320  }
321 
322  //Reset spine height
323  *value_pt = backup;
324 
325  }
326  // Note: Re-assignment of inflow is done on the fly during next loop
327  }
328  }
329 
330  // Final re-assign for the inflow velocities for nodes in this element
331  reassign_inflow();
332 
333  }
334 
335 
336 private:
337 
338 
339 
340  /// \short For all nodes that are located on specified boundary
341  /// re-assign the inflow velocity, using the function pointed to
342  /// by the function pointer
344  {
345  // Loop over all nodes in element -- if they are
346  // on inflow boundary, re-assign their velocities
347  Vector<double> x(2);
348  Vector<double> veloc(2);
349  unsigned n_nod = this->nnode();
350  for (unsigned j=0;j<n_nod;j++)
351  {
352  Node* nod_pt = this->node_pt(j);
353 
354  if(nod_pt->is_on_boundary(Inflow_boundary))
355  {
356 #ifdef PARANOID
357  for (unsigned i=0;i<2;i++)
358  {
359  if (nod_pt->eqn_number(0)>=0)
360  {
361  std::ostringstream error_stream;
362  error_stream
363  << "We're assigning a Dirichlet condition for the "
364  << i << "-th "
365  << "velocity, even though it is not pinned!\n"
366  << "This can't be right! I'm bailing out..."
367  << std::endl;
368 
369  throw OomphLibError(error_stream.str(),
370  OOMPH_CURRENT_FUNCTION,
371  OOMPH_EXCEPTION_LOCATION);
372  }
373  }
374 #endif
375  // Get inflow profile
376  x[0]=nod_pt->x(0);
377  x[1]=nod_pt->x(1);
378  Inflow_fct_pt(x,veloc);
379  nod_pt->set_value(0,veloc[0]);
380  nod_pt->set_value(1,veloc[1]);
381  }
382  }
383  }
384 
385 
386  /// Storage for the external Data that affects the inflow
387  Vector<Data*> Inflow_ext_data;
388 
389  /// \short Storage for the local equation numbers associated the Data
390  /// values that affect the inflow
391  DenseMatrix<int> Inflow_ext_data_eqn;
392 
393  /// Number of the inflow boundary in the global mesh
394  unsigned Inflow_boundary;
395 
396  /// \short Function pointer to the global function that specifies the
397  /// inflow velocity profile on the global mesh boundary Inflow_boundary
398  InflowFctPt Inflow_fct_pt;
399 
400 };
401 
402 
403 
404 // Note: Specialisation must go into namespace!
405 namespace oomph
406 {
407 
408 //=======================================================================
409 /// Face geometry of the Bretherton 2D Crouzeix_Raviart spine elements
410 //=======================================================================
411 template<>
412 class FaceGeometry<BrethertonElement<SpineElement<QCrouzeixRaviartElement<2> > > >: public virtual QElement<1,3>
413 {
414  public:
415  FaceGeometry() : QElement<1,3>() {}
416 };
417 
418 
419 
420 //=======================================================================
421 /// Face geometry of the Bretherton 2D Taylor Hood spine elements
422 //=======================================================================
423 template<>
424 class FaceGeometry<BrethertonElement<SpineElement<QTaylorHoodElement<2> > > >: public virtual QElement<1,3>
425 {
426  public:
427  FaceGeometry() : QElement<1,3>() {}
428 };
429 
430 
431 //=======================================================================
432 /// Face geometry of the Face geometry of the
433 /// the Bretherton 2D Crouzeix_Raviart spine elements
434 //=======================================================================
435 template<>
436 class FaceGeometry<FaceGeometry<BrethertonElement<
437  SpineElement<QCrouzeixRaviartElement<2> > > > >: public virtual PointElement
438 {
439  public:
440  FaceGeometry() : PointElement() {}
441 };
442 
443 
444 
445 //=======================================================================
446 /// Face geometry of face geometry of
447 /// the Bretherton 2D Taylor Hood spine elements
448 //=======================================================================
449 template<>
450 class FaceGeometry<FaceGeometry<BrethertonElement<SpineElement<QTaylorHoodElement<2> > > > >: public virtual PointElement
451 {
452  public:
453  FaceGeometry() : PointElement() {}
454 };
455 
456 
457 }
458 
459 /////////////////////////////////////////////////////////////////////
460 /////////////////////////////////////////////////////////////////////
461 /////////////////////////////////////////////////////////////////////
462 
463 
464 //======================================================================
465 /// Bretherton problem
466 //======================================================================
467 template<class ELEMENT>
468 class BrethertonProblem : public Problem
469 {
470 
471 
472 public:
473 
474  /// Constructor:
476 
477  /// \short Spine heights/lengths are unknowns in the problem so their
478  /// values get corrected during each Newton step. However,
479  /// changing their value does not automatically change the
480  /// nodal positions, so we need to update all of them
482  {
483  // Update
484  Bulk_mesh_pt->node_update();
485 
486  // Apply inflow on boundary 1
487  unsigned ibound=1;
488  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
489 
490  // Coordinate and velocity
491  Vector<double> x(2);
492  Vector<double> veloc(2);
493 
494  // Loop over all nodes
495  for (unsigned inod=0;inod<num_nod;inod++)
496  {
497  // Get nodal position
498  x[0]=mesh_pt()->boundary_node_pt(ibound,inod)->x(0);
499  x[1]=mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
500 
501  // Get inflow profile
503  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,veloc[0]);
504  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,veloc[1]);
505  }
506  }
507 
508  /// \short Update before solve: empty
510 
511  /// \short Update after solve can remain empty, because the update
512  /// is performed automatically after every Newton step.
514 
515  ///Fix pressure value l in element e to value p_value
516  void fix_pressure(const unsigned &e, const unsigned &l,
517  const double &pvalue)
518  {
519  //Fix the pressure at that element
520  dynamic_cast<ELEMENT *>(Bulk_mesh_pt->element_pt(e))->
521  fix_pressure(l,pvalue);
522  }
523 
524 
525  /// \short Activate the dependency of the inflow velocity on the
526  /// spine heights at the outflow
527  void activate_inflow_dependency();
528 
529  /// Run a parameter study; perform specified number of steps
530  void parameter_study(const unsigned& nsteps);
531 
532  /// Doc the solution
533  void doc_solution(DocInfo& doc_info);
534 
535 
536 private:
537 
538  /// Pointer to control element
540 
541  /// Trace file
542  ofstream Trace_file;
543 
544  /// Pointer to bulk mesh
545  BrethertonSpineMesh<ELEMENT,SpineLineFluidInterfaceElement<ELEMENT> >*
547 
548 };
549 
550 
551 //====================================================================
552 /// Problem constructor
553 //====================================================================
554 template<class ELEMENT>
556 {
557 
558  // Number of elements in the deposited film region
559  unsigned nx1=24;
560 
561  // Number of elements in the bottom part of the transition region
562  unsigned nx2=6;
563 
564  // Number of elements in channel region
565  unsigned nx3=12;
566 
567  // Number of elements in the vertical part of the transition region
568  // (=half the number of elements in the liquid filled region ahead
569  // of the finger tip)
570  unsigned nhalf=4;
571 
572  // Number of elements through thickness of deposited film
573  unsigned nh=3;
574 
575  // Thickness of deposited film
576  double h=0.1;
577 
578  // Create wall geom objects
579  GeomObject* lower_wall_pt=new StraightLine(-1.0);
580  GeomObject* upper_wall_pt=new StraightLine( 1.0);
581 
582  // Start coordinate on wall
583  double xi0=-4.0;
584 
585  // End of transition region on wall
586  double xi1=1.5;
587 
588  // End of liquid filled region (inflow) on wall
589  double xi2=5.0;
590 
591  //Now create the mesh
592  Bulk_mesh_pt = new BrethertonSpineMesh<ELEMENT,
593  SpineLineFluidInterfaceElement<ELEMENT> >
594  (nx1,nx2,nx3,nh,nhalf,h,lower_wall_pt,upper_wall_pt,xi0,0.0,xi1,xi2);
595 
596  // Make bulk mesh the global mesh
597  mesh_pt()=Bulk_mesh_pt;
598 
599  // Store the control element
600  Control_element_pt=Bulk_mesh_pt->control_element_pt();
601 
602 
603  // Set pointers to quantities that determine the inflow profile
604 
605  // Film thickness at outflow on lower wall:
607  Bulk_mesh_pt->spine_pt(0)->spine_height_pt()->value_pt(0);
608 
609  // Film thickness at outflow on upper wall:
610  unsigned last_spine=Bulk_mesh_pt->nfree_surface_spines()-1;
612  Bulk_mesh_pt->spine_pt(last_spine)->spine_height_pt()->value_pt(0);
613 
614  // Current y-position on lower wall at inflow
615  unsigned ibound=1;
617  Bulk_mesh_pt->boundary_node_pt(ibound,0)->x_pt(0,1);
618 
619  // Current y-position on upper wall at inflow
620  unsigned nnod=Bulk_mesh_pt->nboundary_node(ibound);
622  Bulk_mesh_pt->boundary_node_pt(ibound,nnod-1)->x_pt(0,1);
623 
624  // Activate dependency of inflow on spine heights at outflow
625  activate_inflow_dependency();
626 
627  // Set the boundary conditions for this problem: All nodes are
628  // free by default -- just pin the ones that have Dirichlet conditions
629  // here
630 
631  // No slip on boundaries 0 1 and 2
632  for(unsigned ibound=0;ibound<=2;ibound++)
633  {
634  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
635  for (unsigned inod=0;inod<num_nod;inod++)
636  {
637  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
638  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
639  }
640  }
641 
642  // Uniform, parallel outflow on boundaries 3 and 5
643  for(unsigned ibound=3;ibound<=5;ibound+=2)
644  {
645  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
646  for (unsigned inod=0;inod<num_nod;inod++)
647  {
648  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
649  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
650  }
651  }
652 
653  // Pin central spine
654  unsigned central_spine=(Bulk_mesh_pt->nfree_surface_spines()-1)/2;
655  Bulk_mesh_pt->spine_pt(central_spine)->spine_height_pt()->pin(0);
656 
657 
658  // No slip in moving frame on boundaries 0 and 2
659  for (unsigned ibound=0;ibound<=2;ibound+=2)
660  {
661  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
662  for (unsigned inod=0;inod<num_nod;inod++)
663  {
664  // Parallel flow
665  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,-1.0);
666  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1, 0.0);
667  }
668  }
669 
670 
671  // Parallel, uniform outflow on boundaries 3 and 5
672  for (unsigned ibound=3;ibound<=5;ibound+=2)
673  {
674  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
675  for (unsigned inod=0;inod<num_nod;inod++)
676  {
677  // Parallel inflow
678  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,-1.0);
679  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1, 0.0);
680  }
681  }
682 
683 
684 
685  //Complete the problem setup to make the elements fully functional
686 
687  //Loop over the elements in the layer
688  unsigned n_bulk=Bulk_mesh_pt->nbulk();
689  for(unsigned i=0;i<n_bulk;i++)
690  {
691  //Cast to a fluid element
692  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->
693  bulk_element_pt(i));
694  //Set the Reynolds number, etc
695  el_pt->re_pt() = &Global_Physical_Variables::Re;
696  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
697  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
698  el_pt->g_pt() = &Global_Physical_Variables::G;
699  }
700 
701  //Loop over 1D interface elements and set capillary number
702  unsigned interface_element_pt_range = Bulk_mesh_pt->ninterface_element();
703  for(unsigned i=0;i<interface_element_pt_range;i++)
704  {
705  //Cast to a interface element
706  SpineLineFluidInterfaceElement<ELEMENT>* el_pt =
707  dynamic_cast<SpineLineFluidInterfaceElement<ELEMENT>*>
708  (Bulk_mesh_pt->interface_element_pt(i));
709 
710  //Set the Capillary number
711  el_pt->ca_pt() = &Global_Physical_Variables::Ca;
712  }
713 
714  //Update the nodal positions, so that the domain is correct
715  //(and therefore the repeated node test passes)
716  Bulk_mesh_pt->node_update();
717 
718  //Do equation numbering
719  cout << "Number of unknowns: " << assign_eqn_numbers() << std::endl;
720 
721 }
722 
723 
724 
725 //========================================================================
726 /// Activate the dependency of the inflow velocity on the
727 /// spine heights at the outflow
728 //========================================================================
729 template<class ELEMENT>
731 {
732 
733  // Spine heights that affect the inflow
734  Vector<Data*> outflow_spines(2);
735 
736  // First spine
737  outflow_spines[0]=Bulk_mesh_pt->spine_pt(0)->spine_height_pt();
738 
739  // Last proper spine
740  unsigned last_spine=Bulk_mesh_pt->nfree_surface_spines()-1;
741  outflow_spines[1]=Bulk_mesh_pt->spine_pt(last_spine)->spine_height_pt();;
742 
743 
744 
745  /// Loop over elements on inflow boundary (1)
746  unsigned ibound=1;
747  unsigned nel=Bulk_mesh_pt->nboundary_element(ibound);
748  for (unsigned e=0;e<nel;e++)
749  {
750  // Get pointer to element
751  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(Bulk_mesh_pt->
752  boundary_element_pt(ibound,e));
753  // Activate depency on inflow
754  el_pt->activate_inflow_dependency_on_external_data(
755  outflow_spines,ibound,&Global_Physical_Variables::inflow);
756  }
757 
758 }
759 
760 
761 
762 //========================================================================
763 /// Doc the solution
764 //========================================================================
765 template<class ELEMENT>
767 {
768 
769  ofstream some_file;
770  char filename[100];
771 
772  // Number of plot points
773  unsigned npts=5;
774 
775  // Control coordinate: At bubble tip
776  Vector<double> s(2);
777  s[0]=1.0;
778  s[1]=1.0;
779 
780  // Last proper spine
781  unsigned last_spine=Bulk_mesh_pt->nfree_surface_spines()-1;
782 
783  // Doc
784  Trace_file << " " << Global_Physical_Variables::Ca;
785  Trace_file << " " << Bulk_mesh_pt->spine_pt(0)->height();
786  Trace_file << " " << Bulk_mesh_pt->spine_pt(last_spine)->height();
787  Trace_file << " " << 1.3375*pow(Global_Physical_Variables::Ca,2.0/3.0);
788  Trace_file << " " << -Control_element_pt->interpolated_p_nst(s)*
790  Trace_file << " " << 1.0+3.8*pow(Global_Physical_Variables::Ca,2.0/3.0);
791  Trace_file << " " << Control_element_pt->interpolated_u_nst(s,0);
792  Trace_file << " " << Control_element_pt->interpolated_u_nst(s,1);
793  Trace_file << std::endl;
794 
795 
796  // Output solution
797  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
798  doc_info.number());
799  some_file.open(filename);
800  Bulk_mesh_pt->output(some_file,npts);
801  some_file.close();
802 
803 
804  // Output boundaries
805  sprintf(filename,"%s/boundaries%i.dat",doc_info.directory().c_str(),
806  doc_info.number());
807  some_file.open(filename);
808  Bulk_mesh_pt->output_boundaries(some_file);
809  some_file.close();
810 
811 }
812 
813 
814 
815 
816 //=============================================================================
817 /// Parameter study
818 //=============================================================================
819 template<class ELEMENT>
821 {
822 
823  // Increase maximum residual
824  Problem::Max_residuals=100.0;
825 
826  // Set output directory
827  DocInfo doc_info;
828  doc_info.set_directory("RESLT");
829  doc_info.number()=0;
830 
831 
832  // Open trace file
833  char filename[100];
834  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
835  Trace_file.open(filename);
836 
837  Trace_file << "VARIABLES=\"Ca\",";
838  Trace_file << "\"h<sub>bottom</sub>\",\"h<sub>too</sub>\",";
839  Trace_file << "\"h<sub>Bretherton</sub>\",\"p<sub>tip</sub>\",";
840  Trace_file << "\"p<sub>tip (Bretherton)</sub>\",\"u<sub>stag</sub>\",";
841  Trace_file << "\"v<sub>stag</sub>\"";
842  Trace_file << "\"<greek>a</greek><sub>bottom</sub>\",";
843  Trace_file << "\"<greek>a</greek><sub>top</sub>\"";
844  Trace_file << std::endl;
845 
846  // Initial scaling factor for Ca reduction
847  double factor=2.0;
848 
849  //Doc initial solution
850  doc_solution(doc_info);
851 
852 //Loop over the steps
853  for (unsigned step=1;step<=nsteps;step++)
854  {
855  cout << std::endl << "STEP " << step << std::endl;
856 
857  // Newton method tends to converge is initial stays bounded below
858  // a certain tolerance: This is cheaper to check than restarting
859  // the Newton method after divergence (we'd also need to back up
860  // the previous solution)
861  double maxres=100.0;
862  while (true)
863  {
864  cout << "Checking max. res for Ca = "
865  << Global_Physical_Variables::Ca << ". Max. residual: ";
866 
867  // Check the maximum residual
868  DoubleVector residuals;
869  actions_before_newton_solve();
870  actions_before_newton_convergence_check();
871  get_residuals(residuals);
872  double max_res=residuals.max();
873  cout << max_res;
874 
875  // Check what to do
876  if (max_res>maxres)
877  {
878  cout << ". Too big!" << std::endl;
879  // Reset the capillary number
880  Global_Physical_Variables::Ca *= factor;
881  // Reduce the factor
882  factor=1.0+(factor-1.0)/1.5;
883  cout << "New reduction factor: " << factor << std::endl;
884  // Reduce the capillary number
885  Global_Physical_Variables::Ca /= factor;
886  // Try again
887  continue;
888  }
889  // Looks promising: Let's proceed to solve
890  else
891  {
892  cout << ". OK" << std::endl << std::endl;
893  // Break out of the Ca adjustment loop
894  break;
895  }
896  }
897 
898  //Solve
899  cout << "Solving for capillary number: "
900  << Global_Physical_Variables::Ca << std::endl;
901  newton_solve();
902 
903  // Doc solution
904  doc_info.number()++;
905  doc_solution(doc_info);
906 
907  // Reduce the capillary number
909  }
910 
911 }
912 
913 
914 
915 //======================================================================
916 /// Driver code for unsteady two-layer fluid problem. If there are
917 /// any command line arguments, we regard this as a validation run
918 /// and perform only a single step.
919 //======================================================================
920 int main(int argc, char* argv[])
921 {
922 
923 
924  // Store command line arguments
925  CommandLineArgs::setup(argc,argv);
926 
927  // Set physical parameters:
928 
929  //Set direction of gravity: Vertically downwards
932 
933  // Womersley number = Reynolds number (St = 1)
936 
937  // The Capillary number
939 
940  // Re/Fr -- a measure of gravity...
942 
943  //Set up the problem
945 
946  // Self test:
947  problem.self_test();
948 
949  // Number of steps:
950  unsigned nstep;
951  if (CommandLineArgs::Argc>1)
952  {
953  // Validation run: Just one step
954  nstep=2;
955  }
956  else
957  {
958  // Full run otherwise
959  nstep=30;
960  }
961 
962  // Run the parameter study: Perform nstep steps
963  problem.parameter_study(nstep);
964 
965 }
966 
double * Y_up_pt
Pointer to y-position at inflow on the upper wall.
Definition: bretherton.cc:76
BrethertonElement()
Constructor: Call the constructor of the underlying element.
Definition: bretherton.cc:167
Vector< Data * > Inflow_ext_data
Storage for the external Data that affects the inflow.
Definition: bretherton.cc:387
BrethertonProblem()
Constructor:
Definition: bretherton.cc:555
void actions_before_newton_convergence_check()
Spine heights/lengths are unknowns in the problem so their values get corrected during each Newton st...
Definition: bretherton.cc:481
double * H_up_pt
Pointer to film thickness at outflow on the upper wall.
Definition: bretherton.cc:70
Vector< double > G(2)
Direction of gravity.
void fix_pressure(const unsigned &e, const unsigned &l, const double &pvalue)
Fix pressure value l in element e to value p_value.
Definition: bretherton.cc:516
void activate_inflow_dependency()
Activate the dependency of the inflow velocity on the spine heights at the outflow.
Definition: bretherton.cc:730
Bretherton problem.
Definition: bretherton.cc:468
double ReSt
Womersley = Reynolds times Strouhal.
Definition: bretherton.cc:55
unsigned Inflow_boundary
Number of the inflow boundary in the global mesh.
Definition: bretherton.cc:394
void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Definition: bretherton.cc:265
void reassign_inflow()
For all nodes that are located on specified boundary re-assign the inflow velocity, using the function pointed to by the function pointer.
Definition: bretherton.cc:343
int main(int argc, char *argv[])
Definition: bretherton.cc:920
void inflow(const Vector< double > &x, Vector< double > &veloc)
Definition: bretherton.cc:80
double Re
Reynolds number.
Definition: bretherton.cc:52
void actions_before_newton_solve()
Update before solve: empty.
Definition: bretherton.cc:509
double * Y_lo_pt
Pointer to y-position at inflow on the lower wall.
Definition: bretherton.cc:73
void actions_after_newton_solve()
Update after solve can remain empty, because the update is performed automatically after every Newton...
Definition: bretherton.cc:513
double ReInvFr
Product of Reynolds and Froude number.
Definition: bretherton.cc:58
ELEMENT * Control_element_pt
Pointer to control element.
Definition: bretherton.cc:539
ofstream Trace_file
Trace file.
Definition: bretherton.cc:542
void activate_inflow_dependency_on_external_data(const Vector< Data * > &inflow_ext_data, const unsigned &inflow_boundary, InflowFctPt inflow_fct_pt)
Activate the dependency of the "inflow" on the external data. Pass the vector of pointers to the exte...
Definition: bretherton.cc:175
void assign_local_eqn_numbers(const bool &store_local_dof_pt)
Definition: bretherton.cc:196
void doc_solution(DocInfo &doc_info)
Doc the solution.
Definition: bretherton.cc:766
DenseMatrix< int > Inflow_ext_data_eqn
Storage for the local equation numbers associated the Data values that affect the inflow...
Definition: bretherton.cc:391
BrethertonSpineMesh< ELEMENT, SpineLineFluidInterfaceElement< ELEMENT > > * Bulk_mesh_pt
Pointer to bulk mesh.
Definition: bretherton.cc:546
void parameter_study(const unsigned &nsteps)
Run a parameter study; perform specified number of steps.
Definition: bretherton.cc:820
InflowFctPt Inflow_fct_pt
Function pointer to the global function that specifies the inflow velocity profile on the global mesh...
Definition: bretherton.cc:398
double Ca
Capillary number.
Definition: bretherton.cc:61
double * H_lo_pt
Pointer to film thickness at outflow on the lower wall.
Definition: bretherton.cc:67