collapsible_channel.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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7//LIC//
8//LIC// This library is free software; you can redistribute it and/or
9//LIC// modify it under the terms of the GNU Lesser General Public
10//LIC// License as published by the Free Software Foundation; either
11//LIC// version 2.1 of the License, or (at your option) any later version.
12//LIC//
13//LIC// This library is distributed in the hope that it will be useful,
14//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16//LIC// Lesser General Public License for more details.
17//LIC//
18//LIC// You should have received a copy of the GNU Lesser General Public
19//LIC// License along with this library; if not, write to the Free Software
20//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21//LIC// 02110-1301 USA.
22//LIC//
23//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24//LIC//
25//LIC//====================================================================
26#include <iostream>
27
28// Generic oomph-lib includes
29#include "generic.h"
30
31// Generic oomph-lib includes
32#include "navier_stokes.h"
33
34//Include the mesh
35#include "meshes/collapsible_channel_mesh.h"
36
37using namespace std;
38
39using namespace oomph;
40
41//==========start_of_BL_Squash =========================================
42/// Namespace to define the mapping [0,1] -> [0,1] that re-distributes
43/// nodal points across the channel width.
44//======================================================================
45namespace BL_Squash
46{
47
48 /// Boundary layer width
49 double Delta=0.1;
50
51 /// Fraction of points in boundary layer
52 double Fract_in_BL=0.5;
53
54 /// Mapping [0,1] -> [0,1] that re-distributes
55 /// nodal points across the channel width
56 double squash_fct(const double& s)
57 {
58 // Default return
59 double y=s;
60 if (s<0.5*Fract_in_BL)
61 {
62 y=Delta*2.0*s/Fract_in_BL;
63 }
64 else if (s>1.0-0.5*Fract_in_BL)
65 {
66 y=2.0*Delta/Fract_in_BL*s+1.0-2.0*Delta/Fract_in_BL;
67 }
68 else
69 {
70 y=(1.0-2.0*Delta)/(1.0-Fract_in_BL)*s+
71 (Delta-0.5*Fract_in_BL)/(1.0-Fract_in_BL);
72 }
73
74 return y;
75 }
76}// end of BL_Squash
77
78
79
80
81
82//===============start_of_oscillating_wall=================================
83/// Straight, horizontal channel wall at \f$ y=H \f$ deforms into an
84/// oscillating parabola. The amplitude of the oscillation
85/// \f$ A \f$ and its period is \f$ T \f$.
86/// The position vector to a point on the wall, parametrised by
87/// the Lagrangian coordinate \f$ \zeta \in [0,L]\f$, is therefore given by
88/// \f[ {\bf R}(\zeta,t) =
89/// \left(
90/// \begin{array}{c}
91/// L_{up} + \zeta \\ 1
92/// \end{array}
93/// \right)
94/// + A
95/// \left(
96/// \begin{array}{l}
97/// - B \sin\left(\frac{2\pi}{L_{collapsible}}\zeta\right) \\ \left(
98/// \frac{2}{L_{collapsible}}\right)^2 \zeta \ (L-\zeta)
99/// \end{array}
100/// \right)
101/// \ \sin\left(\frac{2\pi t}{T}\right)
102/// \ {\cal R}(t)
103/// \f]
104/// The parameter \f$ B \f$ is zero by default. If it is set to a nonzero
105/// value, the material particles on the wall also perform some
106/// horizontal motion. The "ramp" function
107/// \f[
108/// {\cal R}(t) = \left\{
109/// \begin{array}{ll}
110/// \frac{1}{2}(1-\cos(\pi t/T)) & \mbox{for $t<T$} \\ 1 & \mbox{for $t \ge T$}
111/// \end{array}
112/// \right.
113/// \f]
114/// provides a "smooth" startup of the oscillation.
115//=========================================================================
116class OscillatingWall : public GeomObject
117{
118
119public:
120
121 /// Constructor : It's a 2D object, parametrised by
122 /// one Lagrangian coordinate. Arguments: height at ends, x-coordinate of
123 /// left end, length, amplitude of deflection, period of oscillation, and
124 /// pointer to time object
125 OscillatingWall(const double& h, const double& x_left, const double& l,
126 const double& a, const double& period, Time* time_pt) :
127 GeomObject(1,2), H(h), Length(l), X_left(x_left), A(a), B(0.0), T(period),
128 Time_pt(time_pt)
129 {}
130
131 /// Destructor: Empty
132 ~OscillatingWall(){}
133
134/// Access function to the amplitude
135 double& amplitude(){return A;}
136
137/// Access function to the period
138 double& period(){return T;}
139
140 /// Position vector at Lagrangian coordinate zeta
141 /// at time level t.
142 void position(const unsigned& t, const Vector<double>&zeta,
143 Vector<double>& r) const
144 {
145 using namespace MathematicalConstants;
146
147 // Smoothly ramp up the oscillation during the first period
148 double ramp=1.0;
149 if (Time_pt->time(t)<T)
150 {
151 ramp=0.5*(1.0-cos(Pi*Time_pt->time(t)/T));
152 }
153
154 // Position vector
155 r[0] = zeta[0]+X_left
156 -B*A*sin(2.0*3.14159*zeta[0]/Length)*
157 sin(2.0*Pi*(Time_pt->time(t))/T)*ramp;
158
159 r[1] = H+A*((Length-zeta[0])*zeta[0])/pow(0.5*Length,2)*
160 sin(2.0*Pi*(Time_pt->time(t))/T)*ramp;
161
162 } // end of "unsteady" version
163
164
165 /// "Current" position vector at Lagrangian coordinate zeta
166 void position(const Vector<double>&zeta, Vector<double>& r) const
167 {
168 position (0, zeta, r);
169 }
170
171 /// Number of geometric Data in GeomObject: None.
172 unsigned ngeom_data() const {return 0;}
173
174private:
175
176 /// Height at ends
177 double H;
178
179 /// Length
180 double Length;
181
182 /// x-coordinate of left end
183 double X_left;
184
185 /// Amplitude of oscillation
186 double A;
187
188 /// Relative amplitude of horizontal wall motion
189 double B;
190
191 /// Period of the oscillations
192 double T;
193
194 /// Pointer to the global time object
195 Time* Time_pt;
196
197}; // end of oscillating wall
198
199
200
201
202//====start_of_Global_Physical_Variables================
203/// Namespace for phyical parameters
204//======================================================
205namespace Global_Physical_Variables
206{
207 /// Reynolds number
208 double Re=50.0;
209
210 /// Womersley = Reynolds times Strouhal
211 double ReSt=50.0;
212
213 /// Default pressure on the left boundary
214 double P_up=0.0;
215
216 /// Traction required at the left boundary
217 void prescribed_traction(const double& t,
218 const Vector<double>& x,
219 const Vector<double>& n,
220 Vector<double>& traction)
221 {
222 traction.resize(2);
223 traction[0]=P_up;
224 traction[1]=0.0;
225 }
226
227} // end of Global_Physical_Variables
228
229
230
231//=======start_of_problem_class=======================================
232/// Problem class
233//====================================================================
234template <class ELEMENT>
235class CollapsibleChannelProblem : public Problem
236{
237
238 public :
239
240 /// Constructor : the arguments are the number of elements,
241 /// the length of the domain and the amplitude and period of
242 /// the oscillations
243 CollapsibleChannelProblem(const unsigned& nup,
244 const unsigned& ncollapsible,
245 const unsigned& ndown,
246 const unsigned& ny,
247 const double& lup,
248 const double& lcollapsible,
249 const double& ldown,
250 const double& ly,
251 const double& amplitude,
252 const double& period);
253
254 /// Empty destructor
255 ~CollapsibleChannelProblem() {}
256
257 /// Access function for the specific mesh
258 RefineableCollapsibleChannelMesh<ELEMENT>* bulk_mesh_pt()
259 {
260
261 // Upcast from pointer to the Mesh base class to the specific
262 // element type that we're using here.
263 return dynamic_cast<RefineableCollapsibleChannelMesh<ELEMENT>*>
264 (Bulk_mesh_pt);
265
266 } // end of access to bulk mesh
267
268 /// Update the problem specs before solve (empty)
269 void actions_before_newton_solve(){}
270
271 /// Update the problem after solve (empty)
272 void actions_after_newton_solve(){}
273
274 /// Update the velocity boundary condition on the moving wall
275 void actions_before_implicit_timestep();
276
277 /// Actions before adapt: Wipe the mesh of prescribed traction elements
278 void actions_before_adapt();
279
280 /// Actions after adapt: Rebuild the mesh of prescribed traction elements
281 void actions_after_adapt();
282
283 /// Apply initial conditions
284 void set_initial_condition();
285
286 /// Doc the solution
287 void doc_solution(DocInfo& doc_info, ofstream& trace_file);
288
289private :
290
291 /// Create the prescribed traction elements on boundary b
292 /// of the bulk mesh and stick them into the surface mesh.
293 void create_traction_elements(const unsigned &b,
294 Mesh* const &bulk_mesh_pt,
295 Mesh* const &surface_mesh_pt);
296
297 /// Delete prescribed traction elements from the surface mesh
298 void delete_traction_elements(Mesh* const &surface_mesh_pt);
299
300 /// Number of elements in the x direction in the upstream part of the channel
301 unsigned Nup;
302
303 /// Number of elements in the x direction in the "collapsible"
304 /// part of the channel
305 unsigned Ncollapsible;
306
307 /// Number of elements in the x direction in the downstream part of the channel
308 unsigned Ndown;
309
310 /// Number of elements across the channel
311 unsigned Ny;
312
313 /// x-length in the upstream part of the channel
314 double Lup;
315
316 /// x-length in the "collapsible" part of the channel
317 double Lcollapsible;
318
319 /// x-length in the downstream part of the channel
320 double Ldown;
321
322 /// Transverse length
323 double Ly;
324
325 /// Pointer to the geometric object that parametrises the "collapsible" wall
326 OscillatingWall* Wall_pt;
327
328 /// Pointer to the "bulk" mesh
329 RefineableCollapsibleChannelMesh<ELEMENT>* Bulk_mesh_pt;
330
331 /// Pointer to the "surface" mesh that contains the applied traction
332 /// elements
333 Mesh* Surface_mesh_pt;
334
335 /// Pointer to the left control node
336 Node* Left_node_pt;
337
338 /// Pointer to right control node
339 Node* Right_node_pt;
340
341}; // end of problem class
342
343
344
345
346//===start_of_constructor=======================================
347/// Constructor for the collapsible channel problem
348//===============================================================
349template <class ELEMENT>
350CollapsibleChannelProblem<ELEMENT>::CollapsibleChannelProblem(
351 const unsigned& nup,
352 const unsigned& ncollapsible,
353 const unsigned& ndown,
354 const unsigned& ny,
355 const double& lup,
356 const double& lcollapsible,
357 const double& ldown,
358 const double& ly,
359 const double& amplitude,
360 const double& period)
361{
362 // Number of elements
363 Nup=nup;
364 Ncollapsible=ncollapsible;
365 Ndown=ndown;
366 Ny=ny;
367
368 // Lengths of domain
369 Lup=lup;
370 Lcollapsible=lcollapsible;
371 Ldown=ldown;
372 Ly=ly;
373
374 // Overwrite maximum allowed residual to accomodate possibly
375 // poor initial guess for solution
376 Problem::Max_residuals=10000;
377
378 // Allocate the timestepper -- this constructs the Problem's
379 // time object with a sufficient amount of storage to store the
380 // previous timsteps.
381 add_time_stepper_pt(new BDF<2>);
382
383 // Parameters for wall object
384 double height=ly;
385 double length=lcollapsible;
386 double x_left=lup;
387
388 //Create the geometric object that represents the wall
389 Wall_pt=new OscillatingWall(height, x_left, length, amplitude, period,
390 time_pt());
391
392 //Build mesh
393 Bulk_mesh_pt = new RefineableCollapsibleChannelMesh<ELEMENT>(
394 nup, ncollapsible, ndown, ny,
395 lup, lcollapsible, ldown, ly,
396 Wall_pt,
397 time_stepper_pt());
398
399
400 // Enable boundary layer squash function?
401#ifdef USE_BL_SQUASH_FCT
402
403 // Set a non-trivial boundary-layer-squash function...
404 Bulk_mesh_pt->bl_squash_fct_pt() = &BL_Squash::squash_fct;
405
406 // ... and update the nodal positions accordingly
407 Bulk_mesh_pt->node_update();
408
409#endif
410 // end of boundary layer squash function
411
412
413 // Create "surface mesh" that will contain only the prescribed-traction
414 // elements at the inflow. The default constructor just creates the mesh
415 // without giving it any elements, nodes, etc.
416 Surface_mesh_pt = new Mesh;
417
418 // Create prescribed-traction elements from all elements that are
419 // adjacent to boundary 5 (inflow boundary), and add them to the surface mesh.
420 create_traction_elements(5,Bulk_mesh_pt,Surface_mesh_pt);
421
422 // Add the two sub meshes to the problem
423 add_sub_mesh(Bulk_mesh_pt);
424 add_sub_mesh(Surface_mesh_pt);
425
426 // Combine all submeshes added so far into a single Mesh
427 build_global_mesh();
428
429 //Set errror estimator for bulk mesh
430 Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
431 dynamic_cast<RefineableCollapsibleChannelMesh<ELEMENT>*>
432 (Bulk_mesh_pt)->spatial_error_estimator_pt()=error_estimator_pt;
433
434
435 // Loop over the elements to set up element-specific
436 // things that cannot be handled by constructor
437 unsigned n_element=Bulk_mesh_pt->nelement();
438 for(unsigned e=0;e<n_element;e++)
439 {
440 // Upcast from GeneralisedElement to the present element
441 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
442
443 //Set the Reynolds number
444 el_pt->re_pt() = &Global_Physical_Variables::Re;
445
446 // Set the Womersley number
447 el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
448
449 } // end loop over bulk elements
450
451
452
453 // Loop over the traction elements to pass pointer to prescribed
454 // traction function
455 unsigned n_el=Surface_mesh_pt->nelement();
456 for(unsigned e=0;e<n_el;e++)
457 {
458 // Upcast from GeneralisedElement to NavierStokes traction element
459 NavierStokesTractionElement<ELEMENT> *el_pt =
460 dynamic_cast< NavierStokesTractionElement<ELEMENT>*>(
461 Surface_mesh_pt->element_pt(e));
462
463 // Set the pointer to the prescribed traction function
464 el_pt->traction_fct_pt() = &Global_Physical_Variables::prescribed_traction;
465
466 } // end loop over applied traction elements
467
468
469
470
471 //Pin the velocity on the boundaries
472 //x and y-velocities pinned along boundary 0 (bottom boundary) :
473 unsigned ibound=0;
474 unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
475 for (unsigned inod=0;inod<num_nod;inod++)
476 {
477 for(unsigned i=0;i<2;i++)
478 {
479 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
480 }
481 }
482
483
484 //x and y-velocities pinned along boundary 2, 3, 4 (top boundaries) :
485 for(unsigned ib=2;ib<5;ib++)
486 {
487 num_nod= bulk_mesh_pt()->nboundary_node(ib);
488 for (unsigned inod=0;inod<num_nod;inod++)
489 {
490 for(unsigned i=0;i<2;i++)
491 {
492 bulk_mesh_pt()->boundary_node_pt(ib, inod)->pin(i);
493 }
494 }
495 }
496
497 //y-velocity pinned along boundary 1 (right boundary):
498 ibound=1;
499 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
500 for (unsigned inod=0;inod<num_nod;inod++)
501 {
502 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
503 }
504
505
506 //y-velocity pinned along boundary 5 (left boundary):
507 ibound=5;
508 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
509 for (unsigned inod=0;inod<num_nod;inod++)
510 {
511 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
512 }// end of pin_velocity
513
514
515
516 //Select control nodes "half way" up the inflow/outflow cross-sections
517 //--------------------------------------------------------------------
518
519 // Left boundary
520 ibound=5;
521 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
522 unsigned control_nod=num_nod/2;
523 Left_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
524
525 // Right boundary
526 ibound=1;
527 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
528 control_nod=num_nod/2;
529 Right_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
530
531 // Setup equation numbering scheme
532 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
533
534} //end of constructor
535
536
537
538
539//====start_of_doc_solution===================================================
540/// Doc the solution
541//============================================================================
542template <class ELEMENT>
543void CollapsibleChannelProblem<ELEMENT>:: doc_solution(DocInfo& doc_info,
544 ofstream& trace_file)
545{
546 ofstream some_file;
547 char filename[100];
548
549 // Number of plot points
550 unsigned npts;
551 npts=5;
552
553 // Output solution
554 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
555 doc_info.number());
556 some_file.open(filename);
557 bulk_mesh_pt()->output(some_file,npts);
558 some_file.close();
559
560
561 // Get the position of the midpoint on the geometric object
562 Vector<double> zeta(1);
563 zeta[0]=0.5*Lcollapsible;
564 Vector<double> wall_point(2);
565 Wall_pt->position(zeta,wall_point);
566
567 // Write trace file
568 trace_file << time_pt()->time() << " "
569 << wall_point[1] << " "
570 << Left_node_pt->value(0) << " "
571 << Right_node_pt->value(0) << " "
572 << std::endl;
573
574 // Output wall shape
575 sprintf(filename,"%s/wall%i.dat",doc_info.directory().c_str(),
576 doc_info.number());
577 some_file.open(filename);
578 unsigned nplot=100;
579 for (unsigned i=0;i<nplot;i++)
580 {
581 zeta[0]=double(i)/double(nplot-1)*Lcollapsible;
582 Wall_pt->position(zeta,wall_point);
583 some_file << wall_point[0] << " "
584 << wall_point[1] << std::endl;
585 }
586 some_file.close();
587
588} // end_of_doc_solution
589
590
591
592
593//==========start_of_create_traction_elements==================================
594/// Create the traction elements
595//============================================================================
596template <class ELEMENT>
597void CollapsibleChannelProblem<ELEMENT>::create_traction_elements(
598 const unsigned &b, Mesh* const &bulk_mesh_pt, Mesh* const &surface_mesh_pt)
599{
600 // How many bulk elements are adjacent to boundary b?
601 unsigned n_element = bulk_mesh_pt->nboundary_element(b);
602
603 // Loop over the bulk elements adjacent to boundary b
604 for(unsigned e=0;e<n_element;e++)
605 {
606 // Get pointer to the bulk element that is adjacent to boundary b
607 ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>
608 (bulk_mesh_pt->boundary_element_pt(b,e));
609
610 //What is the index of the face of element e that lies along boundary b
611 int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
612
613 // Build the corresponding prescribed-traction element
614 NavierStokesTractionElement<ELEMENT>* traction_element_pt =
615 new NavierStokesTractionElement<ELEMENT>(bulk_elem_pt,face_index);
616
617 //Add the prescribed-flux element to the surface mesh
618 surface_mesh_pt->add_element_pt(traction_element_pt);
619
620 } //end of loop over bulk elements adjacent to boundary b
621
622} // end of create_traction_elements
623
624
625//============start_of_delete_traction_elements==============================
626/// Delete traction elements and wipe the surface mesh
627//=======================================================================
628template<class ELEMENT>
629void CollapsibleChannelProblem<ELEMENT>::
630delete_traction_elements(Mesh* const &surface_mesh_pt)
631{
632 // How many surface elements are in the surface mesh
633 unsigned n_element = surface_mesh_pt->nelement();
634
635 // Loop over the surface elements
636 for(unsigned e=0;e<n_element;e++)
637 {
638 // Kill surface element
639 delete surface_mesh_pt->element_pt(e);
640 }
641
642 // Wipe the mesh
643 surface_mesh_pt->flush_element_and_node_storage();
644
645} // end of delete_traction_elements
646
647
648
649//=======start_of_apply_initial_conditions===================================
650/// Apply initial conditions: Impulsive start from steady Poiseuille flow
651//============================================================================
652template <class ELEMENT>
653void CollapsibleChannelProblem<ELEMENT>::set_initial_condition()
654{
655
656 // Check that timestepper is from the BDF family
657 if (time_stepper_pt()->type()!="BDF")
658 {
659 std::ostringstream error_stream;
660 error_stream
661 << "Timestepper has to be from the BDF family!\n"
662 << "You have specified a timestepper from the "
663 << time_stepper_pt()->type() << " family" << std::endl;
664
665 throw OomphLibError(error_stream.str(),
666 OOMPH_CURRENT_FUNCTION,
667 OOMPH_EXCEPTION_LOCATION);
668 }
669
670 // Update the mesh
671 bulk_mesh_pt()->node_update();
672
673 // Loop over the nodes to set initial guess everywhere
674 unsigned num_nod = bulk_mesh_pt()->nnode();
675 for (unsigned n=0;n<num_nod;n++)
676 {
677 // Get nodal coordinates
678 Vector<double> x(2);
679 x[0]=bulk_mesh_pt()->node_pt(n)->x(0);
680 x[1]=bulk_mesh_pt()->node_pt(n)->x(1);
681
682 // Assign initial condition: Steady Poiseuille flow
683 bulk_mesh_pt()->node_pt(n)->set_value(0,6.0*(x[1]/Ly)*(1.0-(x[1]/Ly)));
684 bulk_mesh_pt()->node_pt(n)->set_value(1,0.0);
685 }
686
687 // Assign initial values for an impulsive start
688 bulk_mesh_pt()->assign_initial_values_impulsive();
689
690
691} // end of set_initial_condition
692
693
694
695//=====start_of_actions_before_implicit_timestep==============================
696/// Execute the actions before timestep: Update the velocity
697/// boundary condition on the moving wall
698//============================================================================
699template <class ELEMENT>
700void CollapsibleChannelProblem<ELEMENT>::actions_before_implicit_timestep()
701{
702 // Update the domain shape
703 bulk_mesh_pt()->node_update();
704
705 // Moving wall: No slip; this implies that the velocity needs
706 // to be updated in response to wall motion
707 unsigned ibound=3;
708 unsigned num_nod=bulk_mesh_pt()->nboundary_node(ibound);
709 for (unsigned inod=0;inod<num_nod;inod++)
710 {
711 // Which node are we dealing with?
712 Node* node_pt=bulk_mesh_pt()->boundary_node_pt(ibound,inod);
713
714 // Apply no slip
715 FSI_functions::apply_no_slip_on_moving_wall(node_pt);
716 }
717} //end of actions_before_implicit_timestep
718
719
720
721
722//=========start_of_actions_before_adapt==================================
723/// Actions before adapt: Wipe the mesh of prescribed traction elements
724//========================================================================
725template<class ELEMENT>
726void CollapsibleChannelProblem<ELEMENT>::actions_before_adapt()
727{
728 // Kill the traction elements and wipe surface mesh
729 delete_traction_elements(Surface_mesh_pt);
730
731 // Rebuild the global mesh.
732 rebuild_global_mesh();
733
734} // end of actions_before_adapt
735
736
737//==========start_of_actions_after_adapt==================================
738/// Actions after adapt: Rebuild the mesh of prescribed traction elements
739//========================================================================
740template<class ELEMENT>
741void CollapsibleChannelProblem<ELEMENT>::actions_after_adapt()
742{
743 // Create prescribed-flux elements from all elements that are
744 // adjacent to boundary 5 and add them to surface mesh
745 create_traction_elements(5,Bulk_mesh_pt,Surface_mesh_pt);
746
747 // Rebuild the global mesh
748 rebuild_global_mesh();
749
750 // Loop over the traction elements to pass pointer to prescribed traction function
751 unsigned n_element=Surface_mesh_pt->nelement();
752 for(unsigned e=0;e<n_element;e++)
753 {
754 // Upcast from GeneralisedElement to NavierStokesTractionElement element
755 NavierStokesTractionElement<ELEMENT> *el_pt =
756 dynamic_cast<NavierStokesTractionElement<ELEMENT>*>(
757 Surface_mesh_pt->element_pt(e));
758
759 // Set the pointer to the prescribed traction function
760 el_pt->traction_fct_pt() = &Global_Physical_Variables::prescribed_traction;
761 }
762} // end of actions_after_adapt
763
764
765
766//=======start_of_driver_code==================================================
767/// Driver code for an unsteady adaptive collapsible channel problem
768/// with prescribed wall motion. Presence of command line arguments
769/// indicates validation run with coarse resolution and small number of
770/// timesteps.
771//=============================================================================
772int main(int argc, char* argv[])
773{
774
775 // Store command line arguments
776 CommandLineArgs::setup(argc,argv);
777
778 // Reduction in resolution for validation run?
779 unsigned coarsening_factor=1;
780 if (CommandLineArgs::Argc>1)
781 {
782 coarsening_factor=4;
783 }
784
785 // Number of elements in the domain
786 unsigned nup=20/coarsening_factor;
787 unsigned ncollapsible=40/coarsening_factor;
788 unsigned ndown=40/coarsening_factor;
789 unsigned ny=16/coarsening_factor;
790
791 // Length of the domain
792 double lup=5.0;
793 double lcollapsible=10.0;
794 double ldown=10.0;
795 double ly=1.0;
796
797 // Initial amplitude of the wall deformation
798 double amplitude=1.0e-2; // ADJUST
799
800 // Period of oscillation
801 double period=0.45;
802
803 // Pressure/applied traction on the left boundary: This is consistent with
804 // steady Poiseuille flow
805 Global_Physical_Variables::P_up=12.0*(lup+lcollapsible+ldown);
806
807
808 //Set output directory
809 DocInfo doc_info;
810 doc_info.set_directory("RESLT");
811
812 // Open a trace file
813 ofstream trace_file;
814 char filename[100];
815 sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
816 trace_file.open(filename);
817
818 // Build the problem with Crouzeix Raviart Elements
819 CollapsibleChannelProblem<RefineableQCrouzeixRaviartElement<2> >
820 problem(nup, ncollapsible, ndown, ny,
821 lup, lcollapsible, ldown, ly,
822 amplitude,period);
823
824
825 // Number of timesteps per period
826 unsigned nsteps_per_period=40;
827
828 // Number of periods
829 unsigned nperiod=3;
830
831 // Number of timesteps (reduced for validation)
832 unsigned nstep=nsteps_per_period*nperiod;
833 if (CommandLineArgs::Argc>1)
834 {
835 nstep=3;
836 }
837
838 //Timestep:
839 double dt=period/double(nsteps_per_period);
840
841 // Start time
842 double t_min=0.0;
843
844 // Initialise timestep and set initial conditions
845 problem.time_pt()->time()=t_min;
846 problem.initialise_dt(dt);
847 problem.set_initial_condition();
848
849 // Output the initial solution
850 problem.doc_solution(doc_info, trace_file);
851
852 // Step number
853 doc_info.number()++;
854
855
856 // Set targets for spatial adaptivity
857 problem.bulk_mesh_pt()->max_permitted_error()=1.0e-3;
858 problem.bulk_mesh_pt()->min_permitted_error()=1.0e-5;
859
860 // Overwrite with reduced targets for validation run to force
861 // some refinement during the first few timesteps
862 if (CommandLineArgs::Argc>1)
863 {
864 problem.bulk_mesh_pt()->max_permitted_error()=1.0e-4;
865 problem.bulk_mesh_pt()->min_permitted_error()=1.0e-6;
866 }
867
868
869 // First timestep: We may re-assign the initial condition
870 // following any mesh adaptation.
871 bool first=true;
872
873 // Max. number of adaptations during first timestep
874 unsigned max_adapt=10;
875
876 // Timestepping loop
877 for (unsigned istep=0;istep<nstep;istep++)
878 {
879 // Solve the problem
880 problem.unsteady_newton_solve(dt, max_adapt, first);
881
882
883 // Outpt the solution
884 problem.doc_solution(doc_info, trace_file);
885
886 // Step number
887 doc_info.number()++;
888
889 // We've done one step: Don't re-assign the initial conditions
890 // and limit the number of adaptive mesh refinements to one
891 // per timestep.
892 first=false;
893 max_adapt=1;
894 }
895
896 trace_file.close();
897
898} //end of driver code
899
900
901
902// {
903// unsigned n=100;
904// for (unsigned i=0;i<n;i++)
905// {
906// double s=double(i)/double(n-1);
907// cout << s << " " << BL_Squash::squash_fct(s) << std::endl;
908// }
909// pause("done");
910// }
Namespace to define the mapping [0,1] -> [0,1] that re-distributes nodal points across the channel wi...
double squash_fct(const double &s)
Mapping [0,1] -> [0,1] that re-distributes nodal points across the channel width.
double Delta
Boundary layer width.
double Fract_in_BL
Fraction of points in boundary layer.
///////////////////////////////////////////////////////////////// ///////////////////////////////////...