multi_domain_ref_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 a Navier--Stokes
27//mesh to an advection diffusion mesh, giving Boussinesq convection
28
29//Oomph-lib headers, we require the generic, advection-diffusion,
30//and navier-stokes elements.
31#include "generic.h"
32#include "advection_diffusion.h"
33#include "navier_stokes.h"
34#include "multi_physics.h"
35
36
37// Both meshes are the 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//======start_of_namespace============================================
47/// Namespace for the physical parameters in the problem
48//====================================================================
50{
51
52 /// Peclet number (identically one from our non-dimensionalisation)
53 double Peclet=1.0;
54
55 /// 1/Prandtl number
56 double Inverse_Prandtl=1.0;
57
58 /// Rayleigh number, set to be greater than
59 /// the threshold for linear instability
60 double Rayleigh = 1800.0;
61
62 /// Gravity vector
63 Vector<double> Direction_of_gravity(2);
64
65} // end_of_namespace
66
67
68/// ///////////////////////////////////////////////////////////////////
69/// ///////////////////////////////////////////////////////////////////
70/// ///////////////////////////////////////////////////////////////////
71
72
73
74//=======start_of_problem_class=======================================
75/// 2D Convection problem on two rectangular domains, discretised
76/// with refineable Navier-Stokes and Advection-Diffusion elements.
77/// The specific type of element is specified via the template parameters.
78//====================================================================
79template<class NST_ELEMENT,class AD_ELEMENT>
80class RefineableConvectionProblem : public Problem
81{
82
83public:
84
85 /// Constructor
87
88 /// Destructor. Empty
90
91 /// Update the problem specs before solve:
93
94 /// Update the problem after solve (empty)
96
97 /// Access function to the NST mesh.
98 /// Casts the pointer to the base Mesh object to
99 /// the actual mesh type.
100 RefineableRectangularQuadMesh<NST_ELEMENT>* nst_mesh_pt()
101 {
102 return dynamic_cast<RefineableRectangularQuadMesh<NST_ELEMENT>*>
103 (Nst_mesh_pt);
104 } // end_of_nst_mesh
105
106 /// Access function to the AD mesh.
107 /// Casts the pointer to the base Mesh object to
108 /// the actual mesh type.
109 RefineableRectangularQuadMesh<AD_ELEMENT>* adv_diff_mesh_pt()
110 {
111 return dynamic_cast<RefineableRectangularQuadMesh<AD_ELEMENT>*>
113 } // end_of_ad_mesh
114
115 /// Actions before adapt:(empty)
117
118 /// Actions after adaptation, reset all sources, then
119 /// re-pin a single pressure degree of freedom
121 {
122 //Unpin all the pressures in NST mesh to avoid pinning two pressures
123 RefineableNavierStokesEquations<2>::
124 unpin_all_pressure_dofs(nst_mesh_pt()->element_pt());
125
126 //Pin the zero-th pressure dof in the zero-th element and set
127 // its value to zero
128 fix_pressure(0,0,0.0);
129
130 // Set external elements for the multi-domain solution.
131 Multi_domain_functions::
132 setup_multi_domain_interactions<NST_ELEMENT,AD_ELEMENT>
133 (this,nst_mesh_pt(),adv_diff_mesh_pt());
134
135 } //end_of_actions_after_adapt
136
137
138 /// Fix pressure in element e at pressure dof pdof and set to pvalue
139 void fix_pressure(const unsigned &e, const unsigned &pdof,
140 const double &pvalue)
141 {
142 //Cast to specific element and fix pressure in NST element
143 if (nst_mesh_pt()->nelement()>0)
144 {
145 dynamic_cast<NST_ELEMENT*>(nst_mesh_pt()->element_pt(e))->
146 fix_pressure(pdof,pvalue);
147 }
148 } // end_of_fix_pressure
149
150
151 /// Set the
152 /// boundary condition on the upper wall to be perturbed slightly
153 /// to force the solution into the symmetry broken state.
155
156 /// Set th
157 /// boundary condition on the upper wall to be unperturbed.
159
160 /// Doc the solution.
161 void doc_solution();
162
163private:
164
165 /// DocInfo object
166 DocInfo Doc_info;
167
168 /// Is the boundary condition imperfect or not
170
171protected:
172
173 /// Navier Stokes mesh
174 RefineableRectangularQuadMesh<NST_ELEMENT>* Nst_mesh_pt;
175
176 /// Advection diffusion mesh
177 RefineableRectangularQuadMesh<AD_ELEMENT>* Adv_diff_mesh_pt;
178
179}; // end of problem class
180
181
182//=======start_of_constructor=============================================
183/// Constructor for adaptive thermal convection problem
184//========================================================================
185template<class NST_ELEMENT,class AD_ELEMENT>
188{
189 // Set output directory
190 Doc_info.set_directory("RESLT");
191
192 // # of elements in x-direction
193 unsigned n_x=9;
194
195 // # of elements in y-direction
196 unsigned n_y=8;
197
198 // Domain length in x-direction
199 double l_x=3.0;
200
201 // Domain length in y-direction
202 double l_y=1.0;
203
204 // Build the meshes
205 Nst_mesh_pt =
206 new RefineableRectangularQuadMesh<NST_ELEMENT>(n_x,n_y,l_x,l_y);
207 Adv_diff_mesh_pt =
208 new RefineableRectangularQuadMesh<AD_ELEMENT>(n_x,n_y,l_x,l_y);
209
210 // Create/set error estimator
211 Nst_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
212 Adv_diff_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
213
214 // Set error targets for adaptive refinement
215 Nst_mesh_pt->max_permitted_error()=0.5e-3;
216 Nst_mesh_pt->min_permitted_error()=0.5e-4;
217 Adv_diff_mesh_pt->max_permitted_error()=0.5e-3;
218 Adv_diff_mesh_pt->min_permitted_error()=0.5e-4;
219
220
221 // Set the boundary conditions for this problem: All nodes are
222 // free by default -- only need to pin the ones that have Dirichlet
223 // conditions here
224
225 //Loop over the boundaries of the NST mesh
226 unsigned num_bound = nst_mesh_pt()->nboundary();
227 for(unsigned ibound=0;ibound<num_bound;ibound++)
228 {
229 //Set the maximum index to be pinned (all values by default)
230 unsigned val_max;
231
232 //Loop over the number of nodes on the boundry
233 unsigned num_nod= nst_mesh_pt()->nboundary_node(ibound);
234 for (unsigned inod=0;inod<num_nod;inod++)
235 {
236
237 //If we are on the side-walls, the v-velocity and temperature
238 //satisfy natural boundary conditions, so we only pin the
239 //first value
240 if((ibound==1) || (ibound==3))
241 {
242 val_max=1;
243 }
244 else
245 {
246 val_max=nst_mesh_pt()->boundary_node_pt(ibound,inod)->nvalue();
247 }
248
249 //Loop over the desired values stored at the nodes and pin
250 for(unsigned j=0;j<val_max;j++)
251 {
252 nst_mesh_pt()->boundary_node_pt(ibound,inod)->pin(j);
253 }
254 }
255 }
256
257 // Pin the zero-th pressure value in the zero-th element and
258 // set its value to zero.
259 fix_pressure(0,0,0.0);
260
261 //Loop over the boundaries of the AD mesh
262 num_bound = adv_diff_mesh_pt()->nboundary();
263 for(unsigned ibound=0;ibound<num_bound;ibound++)
264 {
265 //Set the maximum index to be pinned (all values by default)
266 unsigned val_max;
267
268 //Loop over the number of nodes on the boundry
269 unsigned num_nod= adv_diff_mesh_pt()->nboundary_node(ibound);
270 for (unsigned inod=0;inod<num_nod;inod++)
271 {
272 //If we are on the side-walls, the v-velocity and temperature
273 //satisfy natural boundary conditions, so we don't pin anything
274 // in this mesh
275 if ((ibound==1) || (ibound==3))
276 {
277 val_max=0;
278 }
279 else // pin all values
280 {
281 val_max=adv_diff_mesh_pt()->boundary_node_pt(ibound,inod)->nvalue();
282 //Loop over the desired values stored at the nodes and pin
283 for(unsigned j=0;j<val_max;j++)
284 {
285 adv_diff_mesh_pt()->boundary_node_pt(ibound,inod)->pin(j);
286 }
287 }
288 }
289 } // end of loop over AD mesh boundaries
290
291 // Complete the build of all elements so they are fully functional
292
293 // Loop over the elements to set up element-specific
294 // things that cannot be handled by the (argument-free!) ELEMENT
295 // constructor.
296 unsigned n_nst_element = nst_mesh_pt()->nelement();
297 for(unsigned i=0;i<n_nst_element;i++)
298 {
299 // Upcast from GeneralsedElement to the present element
300 NST_ELEMENT *el_pt = dynamic_cast<NST_ELEMENT*>
301 (nst_mesh_pt()->element_pt(i));
302
303 // Set the Reynolds number (1/Pr in our non-dimensionalisation)
305
306 // Set ReSt (also 1/Pr in our non-dimensionalisation)
308
309 // Set the Rayleigh number
310 el_pt->ra_pt() = &Global_Physical_Variables::Rayleigh;
311
312 //Set Gravity vector
314
315 // We can ignore the external geometric data in the "external"
316 // advection diffusion element when computing the Jacobian matrix
317 // because the interaction does not involve spatial gradients of
318 // the temperature (and also because the mesh isn't moving!)
319 el_pt->ignore_external_geometric_data();
320 }
321
322 unsigned n_ad_element = adv_diff_mesh_pt()->nelement();
323 for(unsigned i=0;i<n_ad_element;i++)
324 {
325 // Upcast from GeneralsedElement to the present element
326 AD_ELEMENT *el_pt = dynamic_cast<AD_ELEMENT*>
327 (adv_diff_mesh_pt()->element_pt(i));
328
329 // Set the Peclet number
330 el_pt->pe_pt() = &Global_Physical_Variables::Peclet;
331
332 // Set the Peclet number multiplied by the Strouhal number
333 el_pt->pe_st_pt() =&Global_Physical_Variables::Peclet;
334
335 // We can ignore the external geometric data in the "external"
336 // Navier Stokes element when computing the Jacobian matrix
337 // because the interaction does not involve spatial gradients of
338 // the velocities (and also because the mesh isn't moving!)
339 el_pt->ignore_external_geometric_data();
340
341 } // end of setup for all AD elements
342
343 // combine the submeshes
344 add_sub_mesh(Nst_mesh_pt);
345 add_sub_mesh(Adv_diff_mesh_pt);
346 build_global_mesh();
347
348 // Set external elements for the multi-domain solution.
349 Multi_domain_functions::
350 setup_multi_domain_interactions<NST_ELEMENT,AD_ELEMENT>
351 (this,nst_mesh_pt(),adv_diff_mesh_pt());
352
353 // Setup equation numbering scheme
354 cout << "Number of equations: " << assign_eqn_numbers() << endl;
355
356} // end of constructor
357
358
359//===================start_actions_before_newton_solve====================
360/// Update the problem specs before solve: (Re-)set boundary conditions
361/// to include an imperfection (or not) depending on the control flag.
362//========================================================================
363template<class NST_ELEMENT,class AD_ELEMENT>
365{
366 // Loop over the boundaries on the NST mesh
367 unsigned num_bound = nst_mesh_pt()->nboundary();
368 for(unsigned ibound=0;ibound<num_bound;ibound++)
369 {
370 // Loop over the nodes on boundary
371 unsigned num_nod=nst_mesh_pt()->nboundary_node(ibound);
372 for(unsigned inod=0;inod<num_nod;inod++)
373 {
374 // Get pointer to node
375 Node* nod_pt=nst_mesh_pt()->boundary_node_pt(ibound,inod);
376
377 //Set the number of velocity components
378 unsigned vel_max=2;
379 //If we are on the side walls we only pin the x-velocity.
380 if((ibound==1) || (ibound==3)) {vel_max = 1;}
381 //Set the pinned velocities to zero
382 for(unsigned j=0;j<vel_max;j++) {nod_pt->set_value(j,0.0);}
383
384 //If we are on the top boundary
385 if(ibound==2)
386 {
387 //Add small velocity imperfection if desired
388 if(Imperfect)
389 {
390 //Read out the x position
391 double x = nod_pt->x(0);
392 //Set a sinusoidal perturbation in the vertical velocity
393 //This perturbation is mass conserving
394 double value = sin(2.0*3.141592654*x/3.0);
395 nod_pt->set_value(1,value);
396 }
397 }
398
399 }
400 }
401
402 // Loop over all the boundaries on the AD mesh
403 num_bound=adv_diff_mesh_pt()->nboundary();
404 for(unsigned ibound=0;ibound<num_bound;ibound++)
405 {
406 // Loop over the nodes on boundary
407 unsigned num_nod=adv_diff_mesh_pt()->nboundary_node(ibound);
408 for(unsigned inod=0;inod<num_nod;inod++)
409 {
410 // Get pointer to node
411 Node* nod_pt=adv_diff_mesh_pt()->boundary_node_pt(ibound,inod);
412
413 //If we are on the top boundary, set the temperature
414 //to -0.5 (cooled)
415 if(ibound==2) {nod_pt->set_value(0,-0.5);}
416
417 //If we are on the bottom boundary, set the temperature
418 //to 0.5 (heated)
419 if(ibound==0) {nod_pt->set_value(0,0.5);}
420 }
421 }
422
423
424} // end of actions before solve
425
426
427
428//====================start_of_doc_solution===============================
429/// Doc the solution
430//========================================================================
431template<class NST_ELEMENT,class AD_ELEMENT>
433{
434 //Declare an output stream and filename
435 ofstream some_file;
436 char filename[100];
437
438 // Number of plot points: npts x npts
439 unsigned npts=5;
440
441 // Output Navier-Stokes solution
442 sprintf(filename,"%s/fluid_soln%i.dat",Doc_info.directory().c_str(),
443 Doc_info.number());
444 some_file.open(filename);
445 nst_mesh_pt()->output(some_file,npts);
446 some_file.close();
447
448 // Output advection diffusion solution
449 sprintf(filename,"%s/temperature_soln%i.dat",Doc_info.directory().c_str(),
450 Doc_info.number());
451 some_file.open(filename);
452 adv_diff_mesh_pt()->output(some_file,npts);
453 some_file.close();
454
455
456 Doc_info.number()++;
457} // end of doc
458
459
460//===============start_of_main========================================
461/// Driver code for 2D Boussinesq convection problem with
462/// adaptivity.
463//====================================================================
464int main()
465{
466
467
468 // Set the direction of gravity
471
472 // Create the problem with 2D nine-node refineable elements.
475 RefineableQCrouzeixRaviartElement<2>,
476 RefineableQAdvectionDiffusionElement<2,3> > ,
478 RefineableQAdvectionDiffusionElement<2,3>,
479 RefineableQCrouzeixRaviartElement<2> >
480 > problem;
481
482 // Apply a perturbation to the upper boundary condition to
483 // force the solution into the symmetry-broken state.
484 problem.enable_imperfection();
485
486 //Solve the problem with (up to) two levels of adaptation
487 problem.newton_solve(2);
488
489 //Document the solution
490 problem.doc_solution();
491
492 // Make the boundary conditions perfect and solve again.
493 // Now the slightly perturbed symmetry broken state computed
494 // above is used as the initial condition and the Newton solver
495 // converges to the symmetry broken solution, even without
496 // the perturbation
497 problem.disable_imperfection();
498 problem.newton_solve(2);
499 problem.doc_solution();
500
501} // end of main
502
503
504
505
506
507
508
509
510
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
void actions_after_newton_solve()
Update the problem after solve (empty)
RefineableRectangularQuadMesh< AD_ELEMENT > * adv_diff_mesh_pt()
Access function to the AD mesh. Casts the pointer to the base Mesh object to the actual mesh type.
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 actions_before_newton_solve()
Update the problem specs before solve:
void disable_imperfection()
Set th boundary condition on the upper wall to be unperturbed.
RefineableRectangularQuadMesh< NST_ELEMENT > * nst_mesh_pt()
Access function to the NST mesh. Casts the pointer to the base Mesh object to the actual mesh type.
void actions_before_adapt()
Actions before adapt:(empty)
void actions_after_adapt()
Actions after adaptation, reset all sources, then re-pin a single pressure degree of freedom.
void enable_imperfection()
Set the boundary condition on the upper wall to be perturbed slightly to force the solution into the ...
RefineableRectangularQuadMesh< NST_ELEMENT > * Nst_mesh_pt
Navier Stokes mesh.
RefineableRectangularQuadMesh< AD_ELEMENT > * Adv_diff_mesh_pt
Advection diffusion mesh.
bool Imperfect
Is the boundary condition imperfect or not.
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
int main()
Driver code for 2D Boussinesq convection problem with adaptivity.
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)