refineable_b_convection.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 a multi-physics problem that couples the Navier--Stokes
27//equations to the advection diffusion equations,
28//giving Boussinesq convection
29
30//Oomph-lib headers, we require the generic, advection-diffusion
31//and navier-stokes elements.
32#include "generic.h"
33#include "advection_diffusion.h"
34#include "navier_stokes.h"
35#include "multi_physics.h"
36
37// The mesh is our standard rectangular quadmesh
38#include "meshes/rectangular_quadmesh.h"
39
40// Use the oomph and std namespaces
41using namespace oomph;
42using namespace std;
43
44
45/// ///////////////////////////////////////////////////////////////////
46/// ///////////////////////////////////////////////////////////////////
47/// ///////////////////////////////////////////////////////////////////
48
49
50//======start_of_namespace============================================
51/// Namespace for the physical parameters in the problem
52//====================================================================
54{
55
56 /// Peclet number (identically one from our non-dimensionalisation)
57 double Peclet=1.0;
58
59 /// 1/Prandtl number
60 double Inverse_Prandtl=1.0;
61
62 /// Rayleigh number, set to be greater than
63 /// the threshold for linear instability
64 double Rayleigh = 1800.0;
65
66 /// Gravity vector
67 Vector<double> Direction_of_gravity(2);
68
69} // end_of_namespace
70
71
72/// ///////////////////////////////////////////////////////////////////
73/// ///////////////////////////////////////////////////////////////////
74/// ///////////////////////////////////////////////////////////////////
75
76
77
78//=======start_of_problem_class=======================================
79/// 2D Convection problem on rectangular domain, discretised
80/// with refineable elements. The specific type
81/// of element is specified via the template parameter.
82//====================================================================
83template<class ELEMENT>
84class RefineableConvectionProblem : public Problem
85{
86
87public:
88
89 /// Constructor
91
92 /// Destructor. Empty
94
95 /// Update the problem specs before solve:
97
98 /// Update the problem after solve (empty)
100
101 /// Overloaded version of the problem's access function to
102 /// the mesh. Recasts the pointer to the base Mesh object to
103 /// the actual mesh type.
104 RectangularQuadMesh<ELEMENT>* mesh_pt()
105 {
106 return dynamic_cast<RectangularQuadMesh<ELEMENT>*>(
107 Problem::mesh_pt());
108 } //end of access function to specic mesh
109
110 /// Actions before adapt:(empty)
112
113 /// Actions after adaptation,
114 /// Re-pin a single pressure degree of freedom
116 {
117 //Unpin all the pressures to avoid pinning two pressures
118 RefineableNavierStokesEquations<2>::
119 unpin_all_pressure_dofs(mesh_pt()->element_pt());
120
121 //Pin the zero-th pressure dof in the zero-th element and set
122 // its value to zero
123 fix_pressure(0,0,0.0);
124 }
125
126 /// Fix pressure in element e at pressure dof pdof and set to pvalue
127 void fix_pressure(const unsigned &e, const unsigned &pdof,
128 const double &pvalue)
129 {
130 //Cast to specific element and fix pressure
131 dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
132 fix_pressure(pdof,pvalue);
133 } // end_of_fix_pressure
134
135 /// Set the
136 /// boundary condition on the upper wall to be perturbed slightly
137 /// to force the solution into the symmetry broken state.
139
140 /// Set the
141 /// boundary condition on the upper wall to be unperturbed.
143
144 /// Doc the solution.
146
147private:
148
149 /// DocInfo object
150 DocInfo Doc_info;
151
152 /// Is the boundary condition imperfect or not
153 bool Imperfect;
154
155}; // end of problem class
156
157
158//=======start_of_constructor=============================================
159/// Constructor for adaptive thermal convection problem
160//========================================================================
161template<class ELEMENT>
163RefineableConvectionProblem() : Imperfect(false)
164{
165 // Set output directory
166 Doc_info.set_directory("RESLT");
167
168 // # of elements in x-direction
169 unsigned n_x=9;
170
171 // # of elements in y-direction
172 unsigned n_y=8;
173
174 // Domain length in x-direction
175 double l_x=3.0;
176
177 // Domain length in y-direction
178 double l_y=1.0;
179
180 // Build the mesh
181 RefineableRectangularQuadMesh<ELEMENT>* cast_mesh_pt =
182 new RefineableRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
183
184 //Set the problem's mesh pointer
185 Problem::mesh_pt() = cast_mesh_pt;
186
187
188 // Create/set error estimator
189 cast_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
190
191 // Set error targets for adaptive refinement
192 cast_mesh_pt->max_permitted_error()=0.5e-3;
193 cast_mesh_pt->min_permitted_error()=0.5e-4;
194
195
196 // Set the boundary conditions for this problem: All nodes are
197 // free by default -- only need to pin the ones that have Dirichlet
198 // conditions here
199
200 //Loop over the boundaries
201 unsigned num_bound = mesh_pt()->nboundary();
202 for(unsigned ibound=0;ibound<num_bound;ibound++)
203 {
204 //Set the maximum index to be pinned (all values by default)
205 unsigned val_max=3;
206 //If we are on the side-walls, the v-velocity and temperature
207 //satisfy natural boundary conditions, so we only pin the
208 //first value
209 if((ibound==1) || (ibound==3)) {val_max=1;}
210
211 //Loop over the number of nodes on the boundry
212 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
213 for (unsigned inod=0;inod<num_nod;inod++)
214 {
215 //Loop over the desired values stored at the nodes and pin
216 for(unsigned j=0;j<val_max;j++)
217 {
218 mesh_pt()->boundary_node_pt(ibound,inod)->pin(j);
219 }
220 }
221 }
222
223 // Pin the zero-th pressure value in the zero-th element and
224 // set its value to zero.
225 fix_pressure(0,0,0.0);
226
227 // Complete the build of all elements so they are fully functional
228
229 // Loop over the elements to set up element-specific
230 // things that cannot be handled by the (argument-free!) ELEMENT
231 // constructor.
232 unsigned n_element = mesh_pt()->nelement();
233 for(unsigned i=0;i<n_element;i++)
234 {
235 // Upcast from GeneralsedElement to the present element
236 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
237
238 // Set the Peclet number
239 el_pt->pe_pt() = &Global_Physical_Variables::Peclet;
240
241 // Set the Peclet Strouhal number
242 el_pt->pe_st_pt() = &Global_Physical_Variables::Peclet;
243
244 // Set the Reynolds number (1/Pr in our non-dimensionalisation)
246
247 // Set ReSt (also 1/Pr in our non-dimensionalisation)
249
250 // Set the Rayleigh number
251 el_pt->ra_pt() = &Global_Physical_Variables::Rayleigh;
252
253 //Set Gravity vector
255 }
256
257 // Setup equation numbering scheme
258 cout <<"Number of equations: " << assign_eqn_numbers() << endl;
259
260} // end of constructor
261
262
263//===================start_actions_before_newton_solve===========================
264/// Update the problem specs before solve: (Re-)set boundary conditions
265/// to include an imperfection (or not) depending on the control flag.
266//========================================================================
267template<class ELEMENT>
269{
270 // Loop over the boundaries
271 unsigned num_bound = mesh_pt()->nboundary();
272 for(unsigned ibound=0;ibound<num_bound;ibound++)
273 {
274 // Loop over the nodes on boundary
275 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
276 for(unsigned inod=0;inod<num_nod;inod++)
277 {
278 // Get pointer to node
279 Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
280
281 //Set the number of velocity components
282 unsigned vel_max=2;
283 //If we are on the side walls we only pin the x-velocity.
284 if((ibound==1) || (ibound==3)) {vel_max = 1;}
285 //Set the pinned velocities to zero
286 for(unsigned j=0;j<vel_max;j++) {nod_pt->set_value(j,0.0);}
287
288 //If we are on the top boundary
289 if(ibound==2)
290 {
291 //Set the temperature to -0.5 (cooled)
292 nod_pt->set_value(2,-0.5);
293 //Add small velocity imperfection if desired
294 if(Imperfect)
295 {
296 //Read out the x position
297 double x = nod_pt->x(0);
298 //Set a sinusoidal perturbation in the vertical velocity
299 //This perturbation is mass conserving
300 double value = sin(2.0*3.141592654*x/3.0);
301 nod_pt->set_value(1,value);
302 }
303 }
304
305 //If we are on the bottom boundary, set the temperature
306 //to 0.5 (heated)
307 if(ibound==0) {nod_pt->set_value(2,0.5);}
308 }
309 }
310
311} // end of actions before solve
312
313
314
315//====================start_of_doc_solution===============================
316/// Doc the solution
317//========================================================================
318template<class ELEMENT>
320{
321 //Declare an output stream and filename
322 ofstream some_file;
323 char filename[100];
324
325 // Number of plot points: npts x npts
326 unsigned npts=5;
327
328 // Output solution
329 //-----------------
330 sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
331 Doc_info.number());
332 some_file.open(filename);
333 mesh_pt()->output(some_file,npts);
334 some_file.close();
335
336 Doc_info.number()++;
337} // end of doc
338
339
340//===============start_of_main========================================
341/// Driver code for 2D Boussinesq convection problem with
342/// adaptivity.
343//====================================================================
344int main()
345{
346
347 // Set the direction of gravity
350
351 // Create the problem with 2D nine-node refineable elements.
354
355 // Apply a perturbation to the upper boundary condition to
356 // force the solution into the symmetry-broken state.
357 problem.enable_imperfection();
358
359 //Solve the problem with (up to) two levels of adaptation
360 problem.newton_solve(2);
361
362 //Document the solution
363 problem.doc_solution();
364
365 // Make the boundary conditions perfect and solve again.
366 // Now the slightly perturbed symmetry broken state computed
367 // above is used as the initial condition and the Newton solver
368 // converges to the symmetry broken solution, even without
369 // the perturbation
370 problem.disable_imperfection();
371 problem.newton_solve(2);
372 problem.doc_solution();
373
374} // end of main
375
376
377
378
379
380
381
382
383
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
void actions_after_newton_solve()
Update the problem after solve (empty)
void fix_pressure(const unsigned &e, const unsigned &pdof, const double &pvalue)
Fix pressure in element e at pressure dof pdof and set to pvalue.
~RefineableConvectionProblem()
Destructor. Empty.
void actions_before_newton_solve()
Update the problem specs before solve:
void disable_imperfection()
Set the boundary condition on the upper wall to be unperturbed.
void doc_solution()
Doc the solution.
void actions_before_adapt()
Actions before adapt:(empty)
void actions_after_adapt()
Actions after adaptation, Re-pin a single pressure degree of freedom.
RefineableConvectionProblem()
Constructor.
void enable_imperfection()
Set the boundary condition on the upper wall to be perturbed slightly to force the solution into the ...
RectangularQuadMesh< ELEMENT > * mesh_pt()
Overloaded version of the problem's access function to the mesh. Recasts the pointer to the base Mesh...
void actions_before_newton_solve()
Update the problem specs before solve:
bool Imperfect
Is the boundary condition imperfect or not.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Namespace for the physical parameters in the problem.
Vector< double > Direction_of_gravity(2)
Gravity vector.
double Rayleigh
Rayleigh number, set to be greater than the threshold for linear instability.
double Inverse_Prandtl
1/Prandtl number
double Peclet
Peclet number (identically one from our non-dimensionalisation)
int main()
Driver code for 2D Boussinesq convection problem with adaptivity.