two_d_adv_diff_SUPG.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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Driver for a simple 2D adv diff problem with SUPG stabilisation
31 
32 //Generic routines
33 #include "generic.h"
34 
35 // The Poisson equations
36 #include "advection_diffusion.h"
37 
38 // The mesh
39 #include "meshes/rectangular_quadmesh.h"
40 
41 using namespace std;
42 
43 using namespace oomph;
44 
45 
46 
47 //======start_of_namespace============================================
48 /// Namespace for global parameters: Unforced problem with
49 /// boundary values corresponding to a steep tanh step profile
50 /// oriented at 45 degrees across the domain.
51 //====================================================================
52 namespace GlobalPhysicalParameters
53 {
54 
55  /// Peclet number
56  double Peclet=200.0;
57 
58  /// Parameter for steepness of step in boundary values
59  double Alpha=50.0;
60 
61  /// Parameter for angle of step in boundary values: 45 degrees
62  double TanPhi=1.0;
63 
64  /// Some "solution" for assignment of boundary values
65  void get_boundary_values(const Vector<double>& x, Vector<double>& u)
66  {
67  u[0]=tanh(1.0-Alpha*(TanPhi*x[0]-x[1]));
68  }
69 
70  /// Zero source function
71  void source_function(const Vector<double>& x_vect, double& source)
72  {
73  source=0.0;
74  }
75 
76  /// Wind
77  void wind_function(const Vector<double>& x, Vector<double>& wind)
78  {
79  wind[0]=sin(6.0*x[1]);
80  wind[1]=cos(6.0*x[0]);
81  }
82 
83 
84 } // end of namespace
85 
86 //////////////////////////////////////////////////////////////////////
87 //////////////////////////////////////////////////////////////////////
88 //////////////////////////////////////////////////////////////////////
89 
90 
91 //====== start_of_problem_class=======================================
92 /// 2D AdvectionDiffusion problem on rectangular domain, discretised
93 /// with refineable 2D QAdvectionDiffusion elements. The specific type
94 /// of element is specified via the template parameter.
95 //====================================================================
96 template<class ELEMENT>
97 class SUPGAdvectionDiffusionProblem : public Problem
98 {
99 
100 public:
101 
102  /// \short Constructor: Pass pointer to source and wind functions, and
103  /// flag to indicate if stabilisation is to be used.
105  AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt source_fct_pt,
106  AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt wind_fct_pt,
107  const bool& use_stabilisation);
108 
109  /// Destructor. Empty
111 
112  /// \short Update the problem specs before solve: Reset boundary conditions
113  /// to the values from the tanh solution and compute stabilisation
114  /// parameter.
115  void actions_before_newton_solve();
116 
117  /// Update the problem after solve (empty)
119 
120  /// \short Doc the solution.
121  void doc_solution();
122 
123  /// \short Overloaded version of the problem's access function to
124  /// the mesh. Recasts the pointer to the base Mesh object to
125  /// the actual mesh type.
126  RectangularQuadMesh<ELEMENT>* mesh_pt()
127  {
128  return dynamic_cast<RectangularQuadMesh<ELEMENT>*>(
129  Problem::mesh_pt());
130  }
131 
132 private:
133 
134  /// DocInfo object
135  DocInfo Doc_info;
136 
137  /// Pointer to source function
138  AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt Source_fct_pt;
139 
140  /// Pointer to wind function
141  AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt Wind_fct_pt;
142 
143  /// Flag to indicate if stabilisation is to be used
145 
146 }; // end of problem class
147 
148 
149 
150 //=====start_of_constructor===============================================
151 /// \short Constructor for AdvectionDiffusion problem: Pass pointer to
152 /// source function and wind functions and flag to indicate
153 /// if stabilisation is to be used.
154 //========================================================================
155 template<class ELEMENT>
157  AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt source_fct_pt,
158  AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt wind_fct_pt,
159  const bool& use_stabilisation)
160  : Source_fct_pt(source_fct_pt), Wind_fct_pt(wind_fct_pt),
161  Use_stabilisation(use_stabilisation)
162 {
163 
164  // Set output directory
165  if (use_stabilisation)
166  {
167  Doc_info.set_directory("RESLT_stabilised");
168  }
169  else
170  {
171  Doc_info.set_directory("RESLT_unstabilised");
172  }
173 
174  // Setup mesh
175 
176  // # of elements in x-direction
177  unsigned n_x=40;
178 
179  // # of elements in y-direction
180  unsigned n_y=40;
181 
182  // Domain length in x-direction
183  double l_x=1.0;
184 
185  // Domain length in y-direction
186  double l_y=2.0;
187 
188  // Build and assign mesh
189  Problem::mesh_pt() =
190  new RectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
191 
192 
193  // Set the boundary conditions for this problem: All nodes are
194  // free by default -- only need to pin the ones that have Dirichlet
195  // conditions here
196  unsigned num_bound = mesh_pt()->nboundary();
197  for(unsigned ibound=0;ibound<num_bound;ibound++)
198  {
199  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
200  for (unsigned inod=0;inod<num_nod;inod++)
201  {
202  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
203  }
204  } // end loop over boundaries
205 
206  // Complete the build of all elements so they are fully functional
207 
208  // Loop over the elements to set up element-specific
209  // things that cannot be handled by the (argument-free!) ELEMENT
210  // constructor: Pass pointer to source function
211  unsigned n_element = mesh_pt()->nelement();
212  for(unsigned i=0;i<n_element;i++)
213  {
214  // Upcast from GeneralsedElement to the present element
215  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
216 
217  //Set the source function pointer
218  el_pt->source_fct_pt() = Source_fct_pt;
219 
220  //Set the wind function pointer
221  el_pt->wind_fct_pt() = Wind_fct_pt;
222 
223  // Set the Peclet number
224  el_pt->pe_pt() = &GlobalPhysicalParameters::Peclet;
225  }
226 
227  // Setup equation numbering scheme
228  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
229 
230 } // end of constructor
231 
232 
233 //========================================start_of_actions_before_newton_solve===
234 /// Update the problem specs before solve: (Re-)set boundary conditions
235 //========================================================================
236 template<class ELEMENT>
238 {
239  // How many boundaries are there?
240  unsigned num_bound = mesh_pt()->nboundary();
241 
242  //Loop over the boundaries
243  for(unsigned ibound=0;ibound<num_bound;ibound++)
244  {
245  // How many nodes are there on this boundary?
246  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
247 
248  // Loop over the nodes on boundary
249  for (unsigned inod=0;inod<num_nod;inod++)
250  {
251  // Get pointer to node
252  Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
253 
254  // Extract nodal coordinates from node:
255  Vector<double> x(2);
256  x[0]=nod_pt->x(0);
257  x[1]=nod_pt->x(1);
258 
259  // Get boundary value
260  Vector<double> u(1);
262 
263  // Assign the value to the one (and only) nodal value at this node
264  nod_pt->set_value(0,u[0]);
265  }
266  }
267 
268  // Now loop over all elements and set the stabilisation parameter
269  unsigned n_element = mesh_pt()->nelement();
270  for(unsigned i=0;i<n_element;i++)
271  {
272  // Upcast from GeneralsedElement to the present element
273  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
274 
275 
276  // Use stabilisation?
277  if (Use_stabilisation)
278  {
279  //Compute stabilisation parameter
280  el_pt->compute_stabilisation_parameter();
281  }
282  else
283  {
284  //Compute stabilisation parameter
285  el_pt->switch_off_stabilisation();
286  }
287 
288  }
289 
290 } // end of actions before solve
291 
292 
293 
294 //===============start_of_doc=============================================
295 /// Doc the solution
296 //========================================================================
297 template<class ELEMENT>
299 {
300 
301  ofstream some_file;
302  char filename[100];
303 
304  // Number of plot points: npts x npts
305  unsigned npts=5;
306 
307  // Output solution
308  //-----------------
309  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
310  Doc_info.number());
311  some_file.open(filename);
312  mesh_pt()->output(some_file,npts);
313  some_file.close();
314 
315 } // end of doc
316 
317 
318 //===== start_of_main=====================================================
319 /// Driver code for 2D AdvectionDiffusion problem
320 //========================================================================
321 int main()
322 {
323 
324  //Set up the problem with stabilisation
325  {
326  bool use_stabilisation=true;
327 
328  // Create the problem with 2D nine-node elements from the
329  // QAdvectionDiffusionElement family. Pass pointer to
330  // source and wind function.
334  use_stabilisation);
335 
336  // Solve the problem
337  problem.newton_solve();
338 
339  //Output the solution
340  problem.doc_solution();
341 
342  }
343 
344 
345 
346  //Set up the problem without stabilisation
347  {
348 
349  bool use_stabilisation=false;
350 
351  // Create the problem with 2D nine-node elements from the
352  // QAdvectionDiffusionElement family. Pass pointer to
353  // source and wind function.
357  use_stabilisation);
358 
359  // Solve the problem
360  problem.newton_solve();
361 
362  //Output the solution
363  problem.doc_solution();
364 
365  }
366 
367 
368 
369 } // end of main
370 
371 
372 
373 
374 
375 
376 
377 
378 
void get_boundary_values(const Vector< double > &x, Vector< double > &u)
Some "solution" for assignment of boundary values.
~SUPGAdvectionDiffusionProblem()
Destructor. Empty.
void actions_after_newton_solve()
Update the problem after solve (empty)
void actions_before_newton_solve()
Update the problem specs before solve: Reset boundary conditions to the values from the tanh solution...
DocInfo Doc_info
DocInfo object.
void wind_function(const Vector< double > &x, Vector< double > &wind)
Wind.
double TanPhi
Parameter for angle of step in boundary values: 45 degrees.
void source_function(const Vector< double > &x_vect, double &source)
Zero source function.
double Peclet
Peclet number.
void doc_solution()
Doc the solution.
int main()
Driver code for 2D AdvectionDiffusion problem.
bool Use_stabilisation
Flag to indicate if stabilisation is to be used.
SUPGAdvectionDiffusionProblem(AdvectionDiffusionEquations< 2 >::AdvectionDiffusionSourceFctPt source_fct_pt, AdvectionDiffusionEquations< 2 >::AdvectionDiffusionWindFctPt wind_fct_pt, const bool &use_stabilisation)
Constructor: Pass pointer to source and wind functions, and flag to indicate if stabilisation is to b...
RectangularQuadMesh< ELEMENT > * mesh_pt()
Overloaded version of the problem's access function to the mesh. Recasts the pointer to the base Mesh...
AdvectionDiffusionEquations< 2 >::AdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function.
double Alpha
Parameter for steepness of step in boundary values.
AdvectionDiffusionEquations< 2 >::AdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function.