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-2023 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
41 using namespace oomph;
42 using 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 //====================================================================
79 template<class NST_ELEMENT,class AD_ELEMENT>
80 class RefineableConvectionProblem : public Problem
81 {
82 
83 public:
84 
85  /// Constructor
87 
88  /// Destructor. Empty
90 
91  /// Update the problem specs before solve:
92  void actions_before_newton_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>*>
112  (Adv_diff_mesh_pt);
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.
154  void enable_imperfection() {Imperfect = true;}
155 
156  /// Set th
157  /// boundary condition on the upper wall to be unperturbed.
158  void disable_imperfection() {Imperfect = false;}
159 
160  /// Doc the solution.
161  void doc_solution();
162 
163 private:
164 
165  /// DocInfo object
166  DocInfo Doc_info;
167 
168  /// Is the boundary condition imperfect or not
169  bool Imperfect;
170 
171 protected:
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 //========================================================================
185 template<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)
307  el_pt->re_st_pt() = &Global_Physical_Variables::Inverse_Prandtl;
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 //========================================================================
363 template<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 //========================================================================
431 template<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 //====================================================================
464 int 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)
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.
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 actions_before_newton_solve()
Update the problem specs before solve:
void disable_imperfection()
Set th boundary condition on the upper wall to be unperturbed.
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.
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 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)