spine_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// Driver for 2-D Channel with changing width, using Taylor Hood
27// and Crouzeix Raviart elements.
28
29// Generic oomph-lib header
30#include "generic.h"
31
32// Navier Stokes headers
33#include "navier_stokes.h"
34
35// The mesh
36#include "meshes/channel_spine_mesh.h"
37
38using namespace std;
39
40using namespace oomph;
41
42//==start_of_namespace===================================================
43/// Namespace for physical parameters
44//=======================================================================
46{
47
48 /// Reynolds number
49 double Re=100;
50
51} // end_of_namespace
52
53
54
55
56
57/// ////////////////////////////////////////////////////////////////////
58/// ////////////////////////////////////////////////////////////////////
59// Sinusoidal wall as geometric object
60/// ////////////////////////////////////////////////////////////////////
61/// ////////////////////////////////////////////////////////////////////
62
63
64
65//=========================================================================
66/// Geometric object representing a sinusoidal wall, parametrised by
67/// \f[ x = \zeta \f]
68/// \f[ y = H + A \sin\left(\frac{2\pi \left(\zeta-\zeta_{\mbox{min}}\right)}
69/// {\zeta_{\mbox{max}}-\zeta_{\mbox{min}}}
70/// \right)\f]
71//=========================================================================
72class SinusoidalWall : public GeomObject
73{
74
75public:
76
77 /// Constructor: Pass height, amplitude, zeta min and zeta max
78 /// (all are pinned by default)
79 SinusoidalWall(const double& height, const double& amplitude,
80 const double& zeta_min, const double& zeta_max)
81 : GeomObject(1,2)
82 {
83 // Make space for the geometric data
84 Geom_data_pt.resize(1);
85
86 // Create data: Four value, no time dependence, free by default
87 Geom_data_pt[0] = new Data(4);
88
89 // I've created the data, I have to clean up
90 Must_clean_up=true;
91
92 // Pin the data
93 for (unsigned i=0;i<4;i++) {Geom_data_pt[0]->pin(i);}
94
95 // Give it a value: Initial height
96 Geom_data_pt[0]->set_value(0,height);
97 Geom_data_pt[0]->set_value(1,amplitude);
98 Geom_data_pt[0]->set_value(2,zeta_min);
99 Geom_data_pt[0]->set_value(3,zeta_max);
100 }
101
102
103 /// Constructor: One item of geometric data, containing
104 /// four values:
105 /// \code
106 /// Geom_data_pt[0]->value(0) = height
107 /// Geom_data_pt[0]->value(1) = amplitude
108 /// Geom_data_pt[0]->value(2) = zeta_min
109 /// Geom_data_pt[0]->value(3) = zeta_max
110 /// \endcode
111 SinusoidalWall(const Vector<Data*>& geom_data_pt) : GeomObject(1,2)
112 {
113#ifdef PARANOID
114 if(geom_data_pt.size()!=1)
115 {
116 std::ostringstream error_stream;
117 error_stream
118 << "Wrong size for geom_data_pt " << geom_data_pt.size() << std::endl;
119 if (geom_data_pt[0]->nvalue()!=1)
120 {
121 error_stream << "Wrong geom_data_pt[0]->nvalue() "
122 << geom_data_pt[0]->nvalue() << std::endl;
123 }
124
125 throw OomphLibError(error_stream.str(),
126 OOMPH_CURRENT_FUNCTION,
127 OOMPH_EXCEPTION_LOCATION);
128 }
129#endif
130 Geom_data_pt.resize(1);
131 Geom_data_pt[0]=geom_data_pt[0];
132
133 // Data has been created externally: Must not clean up
134 Must_clean_up=false;
135 }
136
137
138
139 /// Destructor: Clean up if necessary
140 ~SinusoidalWall()
141 {
142 // Do I need to clean up?
143 if (Must_clean_up)
144 {
145 delete Geom_data_pt[0];
146 Geom_data_pt[0]=0;
147 }
148 }
149
150 /// Access function for horizontal half axis
151 double& height(){return *Geom_data_pt[0]->value_pt(0);}
152
153 /// Access function for vertical half axis
154 double& amplitude(){return *Geom_data_pt[0]->value_pt(1);}
155
156
157 /// Position vector at Lagrangian coordinate zeta
158 void position(const Vector<double>& zeta, Vector<double>& r) const
159 {
160#ifdef PARANOID
161 if (r.size()!=Ndim)
162 {
163 throw OomphLibError("The vector r has the wrong dimension\n",
164 OOMPH_CURRENT_FUNCTION,
165 OOMPH_EXCEPTION_LOCATION);
166 }
167#endif
168 // Get parameters
169 double H = Geom_data_pt[0]->value(0);
170 double A = Geom_data_pt[0]->value(1);
171 double zeta_min = Geom_data_pt[0]->value(2);
172 double zeta_max = Geom_data_pt[0]->value(3);
173 double Lz = zeta_max-zeta_min;
174
175 // Position Vector
176 r[0] = zeta[0];
177 r[1] = H + A*sin(MathematicalConstants::Pi*(zeta[0]-zeta_min)/Lz);
178 }
179
180
181 /// Parametrised position on object: r(zeta). Evaluated at
182 /// previous timestep. t=0: current time; t>0: previous
183 /// timestep.
184 void position(const unsigned& t, const Vector<double>& zeta,
185 Vector<double>& r) const
186 {
187#ifdef PARANOID
188 if (t>Geom_data_pt[0]->time_stepper_pt()->nprev_values())
189 {
190 std::ostringstream error_stream;
191 error_stream
192 << "t > nprev_values() in SpikedLine: " << t << " "
193 << Geom_data_pt[0]->time_stepper_pt()->nprev_values() << std::endl;
194
195 throw OomphLibError(error_stream.str(),
196 OOMPH_CURRENT_FUNCTION,
197 OOMPH_EXCEPTION_LOCATION);
198 }
199#endif
200
201 // Get parameters
202 double H = Geom_data_pt[0]->value(t,0);
203 double A = Geom_data_pt[0]->value(t,1);
204 double zeta_min = Geom_data_pt[0]->value(t,2);
205 double zeta_max = Geom_data_pt[0]->value(t,3);
206 double Lz = zeta_max-zeta_min;
207
208 // Position Vector at time level t
209 r[0] = zeta[0];
210 r[1] = H + A*sin(MathematicalConstants::Pi*(zeta[0]-zeta_min)/Lz);
211 }
212
213
214
215 /// Derivative of position Vector w.r.t. to coordinates:
216 /// \f$ \frac{dR_i}{d \zeta_\alpha}\f$ = drdzeta(alpha,i).
217 /// Evaluated at current time.
218 virtual void dposition(const Vector<double>& zeta,
219 DenseMatrix<double> &drdzeta) const
220 {
221 // Get parametres
222 double A = Geom_data_pt[0]->value(1);
223 double zeta_min = Geom_data_pt[0]->value(2);
224 double zeta_max = Geom_data_pt[0]->value(3);
225 double Lz = zeta_max-zeta_min;
226
227 // Tangent vector
228 drdzeta(0,0)=1.0;
229 drdzeta(0,1)=MathematicalConstants::Pi*A*
230 cos(MathematicalConstants::Pi*(zeta[0]-zeta_min)/Lz)/Lz;
231 }
232
233
234 /// 2nd derivative of position Vector w.r.t. to coordinates:
235 /// \f$ \frac{d^2R_i}{d \zeta_\alpha d \zeta_\beta}\f$ = ddrdzeta(alpha,beta,i).
236 /// Evaluated at current time.
237 virtual void d2position(const Vector<double>& zeta,
238 RankThreeTensor<double> &ddrdzeta) const
239 {
240 // Get parametres
241 double A = Geom_data_pt[0]->value(1);
242 double zeta_min = Geom_data_pt[0]->value(2);
243 double zeta_max = Geom_data_pt[0]->value(3);
244 double Lz = zeta_max-zeta_min;
245
246 // Derivative of tangent vector
247 ddrdzeta(0,0,0)=0.0;
248 ddrdzeta(0,0,1)=-MathematicalConstants::Pi*MathematicalConstants::Pi*
249 A*sin(MathematicalConstants::Pi*(zeta[0]-zeta_min)/Lz)/(Lz*Lz);
250 }
251
252
253 /// Posn Vector and its 1st & 2nd derivatives
254 /// w.r.t. to coordinates:
255 /// \f$ \frac{dR_i}{d \zeta_\alpha}\f$ = drdzeta(alpha,i).
256 /// \f$ \frac{d^2R_i}{d \zeta_\alpha d \zeta_\beta}\f$ = ddrdzeta(alpha,beta,i).
257 /// Evaluated at current time.
258 virtual void d2position(const Vector<double>& zeta, Vector<double>& r,
259 DenseMatrix<double> &drdzeta,
260 RankThreeTensor<double> &ddrdzeta) const
261 {
262 // Get parametres
263 double H = Geom_data_pt[0]->value(0);
264 double A = Geom_data_pt[0]->value(1);
265 double zeta_min = Geom_data_pt[0]->value(2);
266 double zeta_max = Geom_data_pt[0]->value(3);
267 double Lz = zeta_max-zeta_min;
268
269 // Position Vector
270 r[0] = zeta[0];
271 r[1] = H + A*sin(MathematicalConstants::Pi*(zeta[0]-zeta_min)/Lz);
272
273 // Tangent vector
274 drdzeta(0,0)=1.0;
275 drdzeta(0,1)=MathematicalConstants::Pi*A*
276 cos(MathematicalConstants::Pi*(zeta[0]-zeta_min)/Lz)/Lz;
277
278 // Derivative of tangent vector
279 ddrdzeta(0,0,0)=0.0;
280 ddrdzeta(0,0,1)=-MathematicalConstants::Pi*MathematicalConstants::Pi*
281 A*sin(MathematicalConstants::Pi*(zeta[0]-zeta_min)/Lz)/(Lz*Lz);
282 }
283
284 /// How many items of Data does the shape of the object depend on?
285 unsigned ngeom_data() const {return Geom_data_pt.size();}
286
287 /// Return pointer to the j-th Data item that the object's
288 /// shape depends on
289 Data* geom_data_pt(const unsigned& j) {return Geom_data_pt[j];}
290
291
292private:
293
294 /// Vector of pointers to Data items that affects the object's shape
295 Vector<Data*> Geom_data_pt;
296
297 /// Do I need to clean up?
298 bool Must_clean_up;
299
300};
301
302/// ////////////////////////////////////////////////////////////////////
303/// ////////////////////////////////////////////////////////////////////
304/// ////////////////////////////////////////////////////////////////////
305
306
307
308//==start_of_problem_class============================================
309/// Channel flow through a non-uniform channel whose geometry is defined
310/// by a spine mesh.
311//====================================================================
312template<class ELEMENT>
313class ChannelSpineFlowProblem : public Problem
314{
315
316public:
317
318 /// Constructor
320
321 /// Destructor: (empty)
323
324 /// Update the problem specs before solve.
325 /// Set velocity boundary conditions just to be on the safe side...
327 {
328 // Update the mesh
329 mesh_pt()->node_update();
330
331 } // end_of_actions_before_newton_solve
332
333
334 /// Update the after solve (empty)
336
337 /// Doc the solution
338 void doc_solution(DocInfo& doc_info);
339
340private:
341
342 /// Width of channel
343 double Ly;
344
345
346}; // end_of_problem_class
347
348
349
350//==start_of_constructor==================================================
351/// Constructor for ChannelSpineFlow problem
352//========================================================================
353template<class ELEMENT>
355{
356
357 // Setup mesh
358
359 // # of elements in x-direction in left region
360 unsigned Nx0=3;
361 // # of elements in x-direction in centre region
362 unsigned Nx1=12;
363 // # of elements in x-direction in right region
364 unsigned Nx2=8;
365
366 // # of elements in y-direction
367 unsigned Ny=10;
368
369 // Domain length in x-direction in left region
370 double Lx0=0.5;
371 // Domain length in x-direction in centre region
372 double Lx1=0.7;
373 // Domain length in x-direction in right region
374 double Lx2=1.5;
375
376 // Domain length in y-direction
377 Ly=1.0;
378
379 // Build geometric object that represents the sinusoidal bump on
380 // the upper wall:
381
382 // 40% indendentation
383 double amplitude_upper = -0.4*Ly;
384 // Minimum and maximum coordinates of bump
385 double zeta_min=Lx0;
386 double zeta_max=Lx0+Lx1;
387 GeomObject* UpperWall =
388 new SinusoidalWall(Ly,amplitude_upper,zeta_min,zeta_max);
389
390 // Build and assign mesh -- pass pointer to geometric object
391 // that represents the sinusoidal bump on the upper wall
392 Problem::mesh_pt() = new ChannelSpineMesh<ELEMENT>(Nx0,Nx1,Nx2,Ny,
393 Lx0,Lx1,Lx2,Ly,
394 UpperWall);
395
396
397 // Set the boundary conditions for this problem: All nodes are
398 // free by default -- just pin the ones that have Dirichlet conditions
399 // here: All boundaries are Dirichlet boundaries, except on boundary 1
400 unsigned num_bound = mesh_pt()->nboundary();
401 for(unsigned ibound=0;ibound<num_bound;ibound++)
402 {
403 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
404 for (unsigned inod=0;inod<num_nod;inod++)
405 {
406 if (ibound!=1)
407 {
408 // Loop over values (u and v velocities)
409 for (unsigned i=0;i<2;i++)
410 {
411 mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
412 }
413 }
414 else
415 {
416 // Parallel outflow ==> no-slip
417 mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
418 }
419 }
420 } // end loop over boundaries
421
422
423 // No slip on stationary upper and lower walls (boundaries 0 and 2)
424 // and parallel outflow (boundary 1)
425 for (unsigned ibound=0;ibound<num_bound-1;ibound++)
426 {
427 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
428 for (unsigned inod=0;inod<num_nod;inod++)
429 {
430 if (ibound!=1)
431 {
432 for (unsigned i=0;i<2;i++)
433 {
434 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
435 }
436 }
437 else
438 {
439 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
440 }
441 }
442 }
443
444
445 // Setup parabolic inflow along boundary 3:
446 unsigned ibound=3;
447 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
448 for (unsigned inod=0;inod<num_nod;inod++)
449 {
450 double y=mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
451 // Parallel, parabolic inflow
452 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,y*(Ly-y));
453 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
454 }
455
456 // Find number of elements in mesh
457 unsigned n_element = mesh_pt()->nelement();
458
459 // Loop over the elements to set up element-specific
460 // things that cannot be handled by constructor: Pass
461 // pointer to Reynolds number
462 for(unsigned e=0;e<n_element;e++)
463 {
464 // Upcast from GeneralisedElement to the present element
465 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
466 //Set the Reynolds number
467 el_pt->re_pt() = &Global_Physical_Variables::Re;
468 }
469
470 // Setup equation numbering scheme
471 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
472
473} // end_of_constructor
474
475
476
477
478//==start_of_doc_solution=================================================
479/// Doc the solution
480//========================================================================
481template<class ELEMENT>
483{
484
485 ofstream some_file;
486 char filename[100];
487
488 // Number of plot points
489 unsigned npts=5;
490
491 // Output solution
492 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
493 doc_info.number());
494 some_file.open(filename);
495 mesh_pt()->output(some_file,npts);
496 some_file.close();
497
498} // end_of_doc_solution
499
500
501
502//==start_of_main======================================================
503/// Driver for channel flow problem with spine mesh.
504//=====================================================================
505int main()
506{
507
508 // Set output directory
509 DocInfo doc_info;
510 doc_info.set_directory("RESLT");
511 doc_info.number()=0;
512
513 // Solve problem with Taylor Hood elements
514 //---------------------------------------
515 {
516 //Build problem
518 problem;
519
520 // Solve the problem with automatic adaptation
521 problem.newton_solve();
522
523 //Output solution
524 problem.doc_solution(doc_info);
525 // Step number
526 doc_info.number()++;
527
528 } // end of Taylor Hood elements
529
530
531 // Solve problem with Crouzeix Raviart elements
532 //--------------------------------------------
533 {
534 // Build problem
536 problem;
537
538 // Solve the problem with automatic adaptation
539 problem.newton_solve();
540
541 //Output solution
542 problem.doc_solution(doc_info);
543 // Step number
544 doc_info.number()++;
545
546 } // end of Crouzeix Raviart elements
547
548
549
550} // end_of_main
551
552
553
554
555//====================================================================
556/// Create the files to illustrate the sparse node update operation
557//====================================================================
559{
560
561 // Set output directory
562 DocInfo doc_info;
563 doc_info.set_directory("RESLT");
564
565 // Setup mesh
566
567 // # of elements in x-direction in left region
568 unsigned Nx0=1;
569 // # of elements in x-direction in centre region
570 unsigned Nx1=5;
571 // # of elements in x-direction in right region
572 unsigned Nx2=1;
573
574 // # of elements in y-direction
575 unsigned Ny=5;
576
577 // Domain length in x-direction in left region
578 double Lx0=0.5;
579 // Domain length in x-direction in centre region
580 double Lx1=1.0;
581 // Domain length in x-direction in right region
582 double Lx2=0.5;
583
584 // Domain length in y-direction
585 double Ly=1.0;
586
587 double amplitude_upper = -0.4*Ly;
588 double zeta_min=Lx0;
589 double zeta_max=Lx0+Lx1;
590 GeomObject* UpperWall =
591 new SinusoidalWall(Ly,amplitude_upper,zeta_min,zeta_max);
592
593 // Build and assign mesh
594 ChannelSpineMesh<SpineElement<QTaylorHoodElement<2> > >* mesh_pt =
595 new ChannelSpineMesh<SpineElement<QTaylorHoodElement<2> > >(Nx0,Nx1,Nx2,Ny,
596 Lx0,Lx1,Lx2,Ly,
597 UpperWall);
598
599 // Update *all* nodal positions
600 mesh_pt->node_update();
601
602 // Change the amplitude
603 dynamic_cast<SinusoidalWall*>(mesh_pt->wall_pt())->amplitude()=0.5;
604
605 unsigned count=0;
606 ofstream some_file;
607 char filename[100];
608
609 // Number of plot points
610 unsigned npts=5;
611
612 // Output solution
613 sprintf(filename,"%s/mesh_update%i.dat",doc_info.directory().c_str(),
614 count);
615 count++;
616
617 // Loop over spines in centre
618 unsigned n_node = mesh_pt->nnode();
619 for (unsigned inode=0;inode<n_node;inode++)
620 {
621 SpineNode* node_pt=dynamic_cast<SpineNode*>(
622 mesh_pt->node_pt(inode));
623 if (node_pt->node_update_fct_id()==1)
624 {
625 node_pt->node_update();
626 // Output solution
627 some_file.open(filename);
628 sprintf(filename,"%s/mesh_update%i.dat",doc_info.directory().c_str(),
629 count);
630 count++;
631 mesh_pt->output(some_file,npts);
632 some_file.close();
633 }
634 }
635
636}
637
638
639
640
641
642
643
644
645
646
647
648
649
650
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void doc_solution(DocInfo &doc_info)
Doc the solution.
ChannelSpineFlowProblem()
Constructor.
void actions_after_newton_solve()
Update the after solve (empty)
double Ly
Width of channel.
SimpleSpineMesh< ELEMENT > * mesh_pt()
Overload access to mesh.
void actions_before_newton_solve()
Update the problem specs before solve. Set velocity boundary conditions just to be on the safe side....
~ChannelSpineFlowProblem()
Destructor: (empty)
////////////////////////////////////////////////////////////////////
double A
Amplitude of indentation.
double height(const double &x)
Height of domain.
double H
Undeformed height of domain.
void doc_sparse_node_update()
Create the files to illustrate the sparse node update operation.
int main()
Driver for channel flow problem with spine mesh.