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