two_d_adv_diff_adapt2.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 adaptive mesh refinement
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 /// Namespace for physical parameters and boundary conditions
47 /// "sharp" step
48 //====================================================================
49 namespace TanhSolnForAdvectionDiffusion
50 {
51 
52  /// Peclet number
53  double Peclet=200.0;
54 
55  /// Parameter for steepness of step in tanh profile
56  double Alpha;
57 
58  /// Parameter for angle of step in tanh profile
59  double TanPhi;
60 
61  /// Tanh profile for assignment of boundary conditons as a Vector
62  void tanh_profile(const Vector<double>& x, Vector<double>& u)
63  {
64  u[0]=tanh(1.0-Alpha*(TanPhi*x[0]-x[1]));
65  }
66 
67  /// Tanh profile for assignment of boundary conditons as a Vector
68  void tanh_profile(const Vector<double>& x, double& u)
69  {
70  u=tanh(1.0-Alpha*(TanPhi*x[0]-x[1]));
71  }
72 
73  /// Zero source function
74  void source_function(const Vector<double>& x_vect, double& source)
75  {
76  //double x=x_vect[0];
77  //double y=x_vect[1];
78  source = 0.0;
79  }
80 
81  /// Wind: Recirculating cells
82  void wind_function(const Vector<double>& x, Vector<double>& wind)
83  {
84  wind[0]=sin(6.0*x[1]);
85  wind[1]=cos(6.0*x[0]);
86  }
87 
88 } // end of namespace
89 
90 //////////////////////////////////////////////////////////////////////
91 //////////////////////////////////////////////////////////////////////
92 //////////////////////////////////////////////////////////////////////
93 
94 //====== start_of_problem_class=======================================
95 /// 2D AdvectionDiffusion problem on rectangular domain, discretised
96 /// with refineable 2D QAdvectionDiffusion elements. The specific
97 /// type of element is specified via the template parameter.
98 //====================================================================
99 template<class ELEMENT>
100 class RefineableAdvectionDiffusionProblem : public Problem
101 {
102 
103 public:
104 
105  /// Constructor: Pass pointer to source function
107  AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt source_fct_pt,
108  AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt wind_fct_pt);
109 
110  /// Destructor (empty)
112 
113  /// \short Update the problem specs before solve: Reset boundary conditions
114  /// to the values from the exact solution.
115  void actions_before_newton_solve();
116 
117  /// Update the problem after solve (empty)
119 
120  /// \short Doc the solution. DocInfo object stores flags/labels for where the
121  /// output gets written to
122  void doc_solution(DocInfo& doc_info);
123 
124  /// \short Overloaded version of the problem's access function to
125  /// the mesh. Recasts the pointer to the base Mesh object to
126  /// the actual mesh type.
127  RefineableRectangularQuadMesh<ELEMENT>* mesh_pt()
128  {
129  return dynamic_cast<RefineableRectangularQuadMesh<ELEMENT>*>(
130  Problem::mesh_pt());
131  }
132 
133 private:
134 
135  /// Pointer to source function
136  AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt Source_fct_pt;
137 
138  /// Pointer to wind function
139  AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt Wind_fct_pt;
140 
141 }; // end of problem class
142 
143 
144 
145 
146 //=====start_of_constructor===============================================
147 /// Constructor for AdvectionDiffusion problem: Pass pointer to source
148 /// function.
149 //========================================================================
150 template<class ELEMENT>
152  AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt source_fct_pt,
153  AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt wind_fct_pt)
154  : Source_fct_pt(source_fct_pt), Wind_fct_pt(wind_fct_pt)
155 {
156 
157  // Setup mesh
158 
159  // # of elements in x-direction
160  unsigned n_x=4;
161 
162  // # of elements in y-direction
163  unsigned n_y=4;
164 
165  // Domain length in x-direction
166  double l_x=1.0;
167 
168  // Domain length in y-direction
169  double l_y=2.0;
170 
171  // Build and assign mesh
172  Problem::mesh_pt() =
173  new RefineableRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
174 
175  // Create/set error estimator
176  mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;
177 
178  // Set the boundary conditions for this problem: All nodes are
179  // free by default -- only need to pin the ones that have Dirichlet conditions
180  // here.
181  unsigned num_bound = mesh_pt()->nboundary();
182  for(unsigned ibound=0;ibound<num_bound;ibound++)
183  {
184  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
185  for (unsigned inod=0;inod<num_nod;inod++)
186  {
187  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
188  }
189  }
190 
191  // Complete the build of all elements so they are fully functional
192 
193  // Loop over the elements to set up element-specific
194  // things that cannot be handled by the (argument-free!) ELEMENT
195  // constructor: Pass pointer to source function
196  unsigned n_element = mesh_pt()->nelement();
197  for(unsigned i=0;i<n_element;i++)
198  {
199  // Upcast from GeneralsedElement to the present element
200  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
201 
202  //Set the source function pointer
203  el_pt->source_fct_pt() = Source_fct_pt;
204 
205  //Set the wind function pointer
206  el_pt->wind_fct_pt() = Wind_fct_pt;
207 
208  // Set the Peclet number
209  el_pt->pe_pt() = &TanhSolnForAdvectionDiffusion::Peclet;
210  }
211 
212  // Setup equation numbering scheme
213  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
214 
215 } // end of constructor
216 
217 //========================================start_of_actions_before_newton_solve===
218 /// Update the problem specs before solve: (Re-)set boundary conditions
219 /// to the values from the exact solution.
220 //========================================================================
221 template<class ELEMENT>
223 {
224  // How many boundaries are there?
225  unsigned num_bound = mesh_pt()->nboundary();
226 
227  //Loop over the boundaries
228  for(unsigned ibound=0;ibound<num_bound;ibound++)
229  {
230  // How many nodes are there on this boundary?
231  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
232 
233  // Loop over the nodes on boundary
234  for (unsigned inod=0;inod<num_nod;inod++)
235  {
236  // Get pointer to node
237  Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
238 
239  // Extract nodal coordinates from node:
240  Vector<double> x(2);
241  x[0]=nod_pt->x(0);
242  x[1]=nod_pt->x(1);
243 
244  // Compute the value of the exact solution at the nodal point
245  Vector<double> u(1);
247 
248  // Assign the value to the one (and only) nodal value at this node
249  nod_pt->set_value(0,u[0]);
250  }
251  }
252 } // end of actions before solve
253 
254 //===============start_of_doc=============================================
255 /// Doc the solution: doc_info contains labels/output directory etc.
256 //========================================================================
257 template<class ELEMENT>
259  DocInfo& doc_info)
260 {
261 
262  ofstream some_file;
263  char filename[100];
264 
265  // Number of plot points: npts x npts
266  unsigned npts=5;
267 
268  // Output solution
269  //-----------------
270  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
271  doc_info.number());
272  some_file.open(filename);
273  mesh_pt()->output(some_file,npts);
274  some_file.close();
275 
276 } // end of doc
277 
278 //===== start_of_main=====================================================
279 /// Driver code for 2D AdvectionDiffusion problem
280 //========================================================================
281 int main()
282 {
283 
284  //Set up the problem
285  //------------------
286 
287  // Create the problem with 2D nine-node refineable elements from the
288  // RefineableQuadAdvectionDiffusionElement family. Pass pointer to
289  // source and wind function.
290  RefineableAdvectionDiffusionProblem<RefineableQAdvectionDiffusionElement<2,
293 
294  // Create label for output
295  //------------------------
296  DocInfo doc_info;
297 
298  // Set output directory
299  doc_info.set_directory("RESLT");
300 
301  // Step number
302  doc_info.number()=0;
303 
304  // Check if we're ready to go:
305  //----------------------------
306  cout << "\n\n\nProblem self-test ";
307  if (problem.self_test()==0)
308  {
309  cout << "passed: Problem can be solved." << std::endl;
310  }
311  else
312  {
313  throw OomphLibError("Self test failed",
314  OOMPH_CURRENT_FUNCTION,
315  OOMPH_EXCEPTION_LOCATION);
316  }
317 
318  // Set the orientation of the "step" to 45 degrees
320 
321  // Choose a large value for the steepness of the "step"
323 
324  // Doc (default) refinement targets
325  //----------------------------------
326  problem.mesh_pt()->doc_adaptivity_targets(cout);
327 
328 
329 
330  // Solve/doc the problem on the initial, very coarse mesh
331  //-------------------------------------------------------
332 
333  // Solve the problem
334  problem.newton_solve();
335 
336  //Output solution
337  problem.doc_solution(doc_info);
338 
339  //Increment counter for solutions
340  doc_info.number()++; // end of Solve/doc very coarse mesh
341 
342  // Now do (up to) four rounds of fully automatic adapation in response to
343  //-----------------------------------------------------------------------
344  // error estimate
345  //---------------
346  unsigned max_solve=4;
347  for (unsigned isolve=0;isolve<max_solve;isolve++)
348  {
349  // Adapt problem/mesh
350  problem.adapt();
351 
352  // Re-solve the problem if the adaptation has changed anything
353  if ((problem.mesh_pt()->nrefined() !=0)||
354  (problem.mesh_pt()->nunrefined()!=0))
355  {
356  problem.newton_solve();
357  }
358  else
359  {
360  cout << "Mesh wasn't adapted --> we'll stop here" << std::endl;
361  break;
362  }
363 
364  //Output solution
365  problem.doc_solution(doc_info);
366 
367  //Increment counter for solutions
368  doc_info.number()++;
369  } // end (up to) four rounds of fully automatic adapation
370 
371 
372 } //end of main
373 
374 
375 
376 
377 
378 
379 
380 
381 
double TanPhi
Parameter for angle of step.
double Alpha
Parameter for steepness of step.
~RefineableAdvectionDiffusionProblem()
Destructor (empty)
void actions_before_newton_solve()
Update the problem specs before solve: Reset boundary conditions to the values from the tanh solution...
void source_function(const Vector< double > &x_vect, double &source)
Source function required to make the solution above an exact solution.
AdvectionDiffusionEquations< 2 >::AdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function.
void wind_function(const Vector< double > &x, Vector< double > &wind)
Wind.
int main()
Driver code for 2D AdvectionDiffusion problem.
RefineableRectangularQuadMesh< ELEMENT > * mesh_pt()
Overloaded version of the problem's access function to the mesh. Recasts the pointer to the base Mesh...
void tanh_profile(const Vector< double > &x, double &u)
Tanh profile for assignment of boundary conditons as a Vector.
void tanh_profile(const Vector< double > &x, Vector< double > &u)
Tanh profile for assignment of boundary conditons as a Vector.
RefineableAdvectionDiffusionProblem(AdvectionDiffusionEquations< 2 >::AdvectionDiffusionSourceFctPt source_fct_pt, AdvectionDiffusionEquations< 2 >::AdvectionDiffusionWindFctPt wind_fct_pt)
Constructor: Pass pointer to source and wind functions.
void actions_after_newton_solve()
Update the problem after solve (empty)
AdvectionDiffusionEquations< 2 >::AdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function.