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