boussinesq_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//======start_of_namespace============================================
46/// Namespace for the physical parameters in the problem
47//====================================================================
49{
50 /// Peclet number (identically one from our non-dimensionalisation)
51 double Peclet=1.0;
52
53 /// 1/Prandtl number
54 double Inverse_Prandtl=1.0;
55
56 /// Rayleigh number, set to be greater than
57 /// the threshold for linear instability
58 double Rayleigh = 1800.0;
59
60 /// Gravity vector
61 Vector<double> Direction_of_gravity(2);
62
63} // end_of_namespace
64
65/// ///////////////////////////////////////////////////////////////////
66/// ///////////////////////////////////////////////////////////////////
67/// ///////////////////////////////////////////////////////////////////
68
69//====== start_of_problem_class=======================================
70/// 2D Convection problem on rectangular domain, discretised
71/// with refineable elements. The specific type
72/// of element is specified via the template parameter.
73//====================================================================
74template<class ELEMENT>
75class ConvectionProblem : public Problem
76{
77
78public:
79
80 /// Constructor
82
83 /// Destructor. Empty
85
86 /// Update the problem specs before solve (empty)
88
89 /// Update the problem after solve (empty)
91
92 /// Actions before adapt:(empty)
94
95 /// Actions before the timestep (update the the time-dependent
96 /// boundary conditions)
98 {
99 set_boundary_conditions(time_pt()->time());
100 }
101
102 /// Fix pressure in element e at pressure dof pdof and set to pvalue
103 void fix_pressure(const unsigned &e, const unsigned &pdof,
104 const double &pvalue)
105 {
106 //Cast to specific element and fix pressure
107 dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
108 fix_pressure(pdof,pvalue);
109 } // end_of_fix_pressure
110
111 /// Doc the solution.
112 void doc_solution();
113
114 /// Set the boundary conditions
115 void set_boundary_conditions(const double &time);
116
117 /// Overloaded version of the problem's access function to
118 /// the mesh. Recasts the pointer to the base Mesh object to
119 /// the actual mesh type.
120 RectangularQuadMesh<ELEMENT>* mesh_pt()
121 {
122 return dynamic_cast<RectangularQuadMesh<ELEMENT>*>(
123 Problem::mesh_pt());
124 }
125
126private:
127
128 /// DocInfo object
129 DocInfo Doc_info;
130
131}; // end of problem class
132
133//===========start_of_constructor=========================================
134/// Constructor for convection problem
135//========================================================================
136template<class ELEMENT>
138{
139 //Allocate a timestepper
140 add_time_stepper_pt(new BDF<2>);
141
142 // Set output directory
143 Doc_info.set_directory("RESLT");
144
145 // # of elements in x-direction
146 unsigned n_x=8;
147
148 // # of elements in y-direction
149 unsigned n_y=8;
150
151 // Domain length in x-direction
152 double l_x=3.0;
153
154 // Domain length in y-direction
155 double l_y=1.0;
156
157 // Build a standard rectangular quadmesh
158 Problem::mesh_pt() =
159 new RectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y,time_stepper_pt());
160
161 // Set the boundary conditions for this problem: All nodes are
162 // free by default -- only need to pin the ones that have Dirichlet
163 // conditions here
164
165 //Loop over the boundaries
166 unsigned num_bound = mesh_pt()->nboundary();
167 for(unsigned ibound=0;ibound<num_bound;ibound++)
168 {
169 //Set the maximum index to be pinned (all values by default)
170 unsigned val_max=3;
171 //If we are on the side-walls, the v-velocity and temperature
172 //satisfy natural boundary conditions, so we only pin the
173 //first value
174 if((ibound==1) || (ibound==3)) {val_max=1;}
175
176 //Loop over the number of nodes on the boundry
177 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
178 for (unsigned inod=0;inod<num_nod;inod++)
179 {
180 //Loop over the desired values stored at the nodes and pin
181 for(unsigned j=0;j<val_max;j++)
182 {
183 mesh_pt()->boundary_node_pt(ibound,inod)->pin(j);
184 }
185 }
186 }
187
188 //Pin the zero-th pressure dof in element 0 and set its value to
189 //zero:
190 fix_pressure(0,0,0.0);
191
192 // Complete the build of all elements so they are fully functional
193
194 // Loop over the elements to set up element-specific
195 // things that cannot be handled by the (argument-free!) ELEMENT
196 // constructor.
197 unsigned n_element = mesh_pt()->nelement();
198 for(unsigned i=0;i<n_element;i++)
199 {
200 // Upcast from GeneralsedElement to the present element
201 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
202
203 // Set the Peclet number
204 el_pt->pe_pt() = &Global_Physical_Variables::Peclet;
205
206 // Set the Peclet number multiplied by the Strouhal number
207 el_pt->pe_st_pt() =&Global_Physical_Variables::Peclet;
208
209 // Set the Reynolds number (1/Pr in our non-dimensionalisation)
211
212 // Set ReSt (also 1/Pr in our non-dimensionalisation)
214
215 // Set the Rayleigh number
216 el_pt->ra_pt() = &Global_Physical_Variables::Rayleigh;
217
218 //Set Gravity vector
220
221 //The mesh is fixed, so we can disable ALE
222 el_pt->disable_ALE();
223 }
224
225 // Setup equation numbering scheme
226 cout <<"Number of equations: " << assign_eqn_numbers() << endl;
227
228} // end of constructor
229
230
231
232//===========start_of_set_boundary_conditions================
233/// Set the boundary conditions as a function of continuous
234/// time
235//===========================================================
236template<class ELEMENT>
238 const double &time)
239{
240 // Loop over the boundaries
241 unsigned num_bound = mesh_pt()->nboundary();
242 for(unsigned ibound=0;ibound<num_bound;ibound++)
243 {
244 // Loop over the nodes on boundary
245 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
246 for(unsigned inod=0;inod<num_nod;inod++)
247 {
248 // Get pointer to node
249 Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
250
251 //Set the number of velocity components
252 unsigned vel_max=2;
253
254 //If we are on the side walls we only set the x-velocity.
255 if((ibound==1) || (ibound==3)) {vel_max = 1;}
256
257 //Set the pinned velocities to zero
258 for(unsigned j=0;j<vel_max;j++) {nod_pt->set_value(j,0.0);}
259
260 //If we are on the top boundary
261 if(ibound==2)
262 {
263 //Set the temperature to -0.5 (cooled)
264 nod_pt->set_value(2,-0.5);
265
266 //Add small velocity imperfection if desired
267 double epsilon = 0.01;
268
269 //Read out the x position
270 double x = nod_pt->x(0);
271
272 //Set a sinusoidal perturbation in the vertical velocity
273 //This perturbation is mass conserving
274 double value = sin(2.0*MathematicalConstants::Pi*x/3.0)*
275 epsilon*time*exp(-time);
276 nod_pt->set_value(1,value);
277 }
278
279 //If we are on the bottom boundary, set the temperature
280 //to 0.5 (heated)
281 if(ibound==0) {nod_pt->set_value(2,0.5);}
282 }
283 }
284} // end_of_set_boundary_conditions
285
286//===============start_doc_solution=======================================
287/// Doc the solution
288//========================================================================
289template<class ELEMENT>
291{
292 //Declare an output stream and filename
293 ofstream some_file;
294 char filename[100];
295
296 // Number of plot points: npts x npts
297 unsigned npts=5;
298
299 // Output solution
300 //-----------------
301 sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
302 Doc_info.number());
303 some_file.open(filename);
304 mesh_pt()->output(some_file,npts);
305 some_file.close();
306
307 Doc_info.number()++;
308} // end of doc
309
310
311//=======start_of_main================================================
312/// Driver code for 2D Boussinesq convection problem
313//====================================================================
314int main(int argc, char **argv)
315{
316
317 // Set the direction of gravity
320
321 //Construct our problem
323
324 // Apply the boundary condition at time zero
325 problem.set_boundary_conditions(0.0);
326
327 //Perform a single steady Newton solve
328 problem.steady_newton_solve();
329
330 //Document the solution
331 problem.doc_solution();
332
333 //Set the timestep
334 double dt = 0.1;
335
336 //Initialise the value of the timestep and set an impulsive start
337 problem.assign_initial_values_impulsive(dt);
338
339 //Set the number of timesteps to our default value
340 unsigned n_steps = 200;
341
342 //If we have a command line argument, perform fewer steps
343 //(used for self-test runs)
344 if(argc > 1) {n_steps = 5;}
345
346 //Perform n_steps timesteps
347 for(unsigned i=0;i<n_steps;++i)
348 {
349 problem.unsteady_newton_solve(dt);
350 problem.doc_solution();
351 }
352
353} // end of main
354
355
356
357
358
359
360
361
362
int main(int argc, char **argv)
Driver code for 2D Boussinesq convection problem.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
ConvectionProblem()
Constructor.
void actions_before_newton_solve()
Update the problem specs before solve (empty)
void actions_before_implicit_timestep()
Actions before the timestep (update the the time-dependent boundary conditions)
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.
void set_boundary_conditions(const double &time)
Set the boundary conditions.
void actions_after_newton_solve()
Update the problem after solve (empty)
~ConvectionProblem()
Destructor. Empty.
void doc_solution()
Doc the solution.
void actions_before_adapt()
Actions before adapt:(empty)
RectangularQuadMesh< ELEMENT > * mesh_pt()
Overloaded version of the problem's access function to the mesh. Recasts the pointer to the base Mesh...
DocInfo Doc_info
DocInfo object.
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)