two_d_unsteady_heat.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 2D unsteady heat problem
27
28//Generic routines
29#include "generic.h"
30
31// The unsteady heat equations
32#include "unsteady_heat.h"
33
34// Mesh
35#include "meshes/rectangular_quadmesh.h"
36
37using namespace std;
38
39using namespace oomph;
40
41using namespace MathematicalConstants;
42
43
44/// ////////////////////////////////////////////////////////////////////
45/// ////////////////////////////////////////////////////////////////////
46/// ////////////////////////////////////////////////////////////////////
47
48//======start_of_ExactSolnForUnsteadyHeat=====================
49/// Namespace for unforced exact solution for UnsteadyHeat equation
50//====================================================================
52{
53
54 /// Decay factor
55 double K=10;
56
57 /// Angle of bump
58 double Phi=1.0;
59
60 /// Exact solution as a Vector
61 void get_exact_u(const double& time, const Vector<double>& x,
62 Vector<double>& u)
63 {
64 double zeta=cos(Phi)*x[0]+sin(Phi)*x[1];
65 u[0]=exp(-K*time)*sin(zeta*sqrt(K));
66 }
67
68 /// Exact solution as a scalar
69 void get_exact_u(const double& time, const Vector<double>& x, double& u)
70 {
71 double zeta=cos(Phi)*x[0]+sin(Phi)*x[1];
72 u=exp(-K*time)*sin(zeta*sqrt(K));
73 }
74
75 /// Source function to make it an exact solution
76 void get_source(const double& time, const Vector<double>& x, double& source)
77 {
78 source = 0.0;
79 }
80
81} // end of ExactSolnForUnsteadyHeat
82
83/// /////////////////////////////////////////////////////////////////////
84/// /////////////////////////////////////////////////////////////////////
85/// /////////////////////////////////////////////////////////////////////
86
87//=====start_of_problem_class=========================================
88/// UnsteadyHeat problem
89//====================================================================
90template<class ELEMENT>
91class UnsteadyHeatProblem : public Problem
92{
93
94public:
95
96 /// Constructor
97 UnsteadyHeatProblem(UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt
98 source_fct_pt);
99
100 /// Destructor (empty)
102
103 /// Update the problem specs after solve (empty)
105
106 /// Update the problem specs before solve (empty)
108
109 /// Update the problem specs after solve (empty)
111
112 /// Update the problem specs before next timestep:
113 /// Set Dirchlet boundary conditions from exact solution.
115
116 /// Set initial condition (incl previous timesteps) according
117 /// to specified function.
119
120 /// Doc the solution
121 void doc_solution(DocInfo& doc_info, ofstream& trace_file);
122
123private:
124
125 /// Pointer to source function
126 UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt Source_fct_pt;
127
128 /// Pointer to control node at which the solution is documented
130
131}; // end of problem class
132
133
134//========start_of_constructor============================================
135/// Constructor for UnsteadyHeat problem in square domain
136//========================================================================
137template<class ELEMENT>
139 UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt source_fct_pt) :
140 Source_fct_pt(source_fct_pt)
141{
142
143 // Allocate the timestepper -- this constructs the Problem's
144 // time object with a sufficient amount of storage to store the
145 // previous timsteps.
146 add_time_stepper_pt(new BDF<2>);
147
148 // Setup parameters for exact solution
149 // -----------------------------------
150
151 // Decay parameter
153
154 // Setup mesh
155 //-----------
156
157 // Number of elements in x and y directions
158 unsigned nx=5;
159 unsigned ny=5;
160
161 // Lengths in x and y directions
162 double lx=1.0;
163 double ly=1.0;
164
165 // Build mesh
166 mesh_pt() = new RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,time_stepper_pt());
167
168 // Choose a control node at which the solution is documented
169 //----------------------------------------------------------
170 // Total number of elements
171 unsigned n_el=mesh_pt()->nelement();
172
173 // Choose an element in the middle
174 unsigned control_el=unsigned(n_el/2);
175
176 // Choose its first node as the control node
177 Control_node_pt=mesh_pt()->finite_element_pt(control_el)->node_pt(0);
178
179 cout << "Recording trace of the solution at: "
180 << Control_node_pt->x(0) << " "
181 << Control_node_pt->x(1) << std::endl;
182
183
184 // Set the boundary conditions for this problem:
185 // ---------------------------------------------
186 // All nodes are free by default -- just pin the ones that have
187 // Dirichlet conditions here.
188 unsigned n_bound = mesh_pt()->nboundary();
189 for(unsigned b=0;b<n_bound;b++)
190 {
191 unsigned n_node = mesh_pt()->nboundary_node(b);
192 for (unsigned n=0;n<n_node;n++)
193 {
194 mesh_pt()->boundary_node_pt(b,n)->pin(0);
195 }
196 } // end of set boundary conditions
197
198
199 // Complete the build of all elements so they are fully functional
200 //----------------------------------------------------------------
201
202 // Find number of elements in mesh
203 unsigned n_element = mesh_pt()->nelement();
204
205 // Loop over the elements to set up element-specific
206 // things that cannot be handled by constructor
207 for(unsigned i=0;i<n_element;i++)
208 {
209 // Upcast from FiniteElement to the present element
210 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
211
212 //Set the source function pointer
213 el_pt->source_fct_pt() = Source_fct_pt;
214 }
215
216 // Do equation numbering
217 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
218
219} // end of constructor
220
221
222
223//=========start of actions_before_implicit_timestep===============================
224/// Actions before timestep: update the domain, then reset the
225/// boundary conditions for the current time.
226//========================================================================
227template<class ELEMENT>
229{
230 // Get current time
231 double time=time_pt()->time();
232
233 //Loop over the boundaries
234 unsigned num_bound = mesh_pt()->nboundary();
235 for(unsigned ibound=0;ibound<num_bound;ibound++)
236 {
237 // Loop over the nodes on boundary
238 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
239 for (unsigned inod=0;inod<num_nod;inod++)
240 {
241 Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
242 double u;
243 Vector<double> x(2);
244 x[0]=nod_pt->x(0);
245 x[1]=nod_pt->x(1);
246 // Get current values of the boundary conditions from the
247 // exact solution
249 nod_pt->set_value(0,u);
250 }
251 }
252} // end of actions_before_implicit_timestep
253
254
255
256//======================start_of_set_initial_condition====================
257/// Set initial condition: Assign previous and current values
258/// from exact solution.
259//========================================================================
260template<class ELEMENT>
262{
263 // Backup time in global Time object
264 double backed_up_time=time_pt()->time();
265
266 // Past history needs to be established for t=time0-deltat, ...
267 // Then provide current values (at t=time0) which will also form
268 // the initial guess for the first solve at t=time0+deltat
269
270 // Vector of exact solution value
271 Vector<double> soln(1);
272 Vector<double> x(2);
273
274 //Find number of nodes in mesh
275 unsigned num_nod = mesh_pt()->nnode();
276
277 // Set continuous times at previous timesteps:
278 // How many previous timesteps does the timestepper use?
279 int nprev_steps=time_stepper_pt()->nprev_values();
280 Vector<double> prev_time(nprev_steps+1);
281 for (int t=nprev_steps;t>=0;t--)
282 {
283 prev_time[t]=time_pt()->time(unsigned(t));
284 }
285
286 // Loop over current & previous timesteps
287 for (int t=nprev_steps;t>=0;t--)
288 {
289 // Continuous time
290 double time=prev_time[t];
291 cout << "setting IC at time =" << time << std::endl;
292
293 // Loop over the nodes to set initial guess everywhere
294 for (unsigned n=0;n<num_nod;n++)
295 {
296 // Get nodal coordinates
297 x[0]=mesh_pt()->node_pt(n)->x(0);
298 x[1]=mesh_pt()->node_pt(n)->x(1);
299
300 // Get exact solution at previous time
302
303 // Assign solution
304 mesh_pt()->node_pt(n)->set_value(t,0,soln[0]);
305
306 // Loop over coordinate directions: Mesh doesn't move, so
307 // previous position = present position
308 for (unsigned i=0;i<2;i++)
309 {
310 mesh_pt()->node_pt(n)->x(t,i)=x[i];
311 }
312 }
313 }
314
315 // Reset backed up time for global timestepper
316 time_pt()->time()=backed_up_time;
317
318} // end of set_initial_condition
319
320
321
322//=======start_of_doc_solution============================================
323/// Doc the solution
324//========================================================================
325template<class ELEMENT>
327doc_solution(DocInfo& doc_info,ofstream& trace_file)
328{
329 ofstream some_file;
330 char filename[100];
331
332 // Number of plot points
333 unsigned npts;
334 npts=5;
335
336
337 cout << std::endl;
338 cout << "=================================================" << std::endl;
339 cout << "Docing solution for t=" << time_pt()->time() << std::endl;
340 cout << "=================================================" << std::endl;
341
342
343 // Output solution
344 //-----------------
345 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
346 doc_info.number());
347 some_file.open(filename);
348 mesh_pt()->output(some_file,npts);
349
350 // Write file as a tecplot text object
351 some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
352 << time_pt()->time() << "\"";
353 // ...and draw a horizontal line whose length is proportional
354 // to the elapsed time
355 some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
356 some_file << "1" << std::endl;
357 some_file << "2" << std::endl;
358 some_file << " 0 0" << std::endl;
359 some_file << time_pt()->time()*20.0 << " 0" << std::endl;
360 some_file.close();
361
362
363 // Output exact solution
364 //----------------------
365 sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
366 doc_info.number());
367 some_file.open(filename);
368 mesh_pt()->output_fct(some_file,npts,time_pt()->time(),
370 some_file.close();
371
372 // Doc error
373 //----------
374 double error,norm;
375 sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
376 doc_info.number());
377 some_file.open(filename);
378 mesh_pt()->compute_error(some_file,
380 time_pt()->time(),
381 error,norm);
382 some_file.close();
383
384 // Doc solution and error
385 //-----------------------
386 cout << "error: " << error << std::endl;
387 cout << "norm : " << norm << std::endl << std::endl;
388
389 // Get exact solution at control node
390 Vector<double> x_ctrl(2);
391 x_ctrl[0]=Control_node_pt->x(0);
392 x_ctrl[1]=Control_node_pt->x(1);
393 double u_exact;
394 ExactSolnForUnsteadyHeat::get_exact_u(time_pt()->time(),x_ctrl,u_exact);
395 trace_file << time_pt()->time() << " "
396 << Control_node_pt->value(0) << " "
397 << u_exact << " "
398 << error << " "
399 << norm << std::endl;
400
401} // end of doc_solution
402
403
404
405/// /////////////////////////////////////////////////////////////////////
406/// /////////////////////////////////////////////////////////////////////
407/// /////////////////////////////////////////////////////////////////////
408
409
410
411//=======start_of_main====================================================
412/// Driver code for unsteady heat equation
413//========================================================================
414int main()
415{
416
417 // Build problem
420
421 // Setup labels for output
422 DocInfo doc_info;
423
424 // Output directory
425 doc_info.set_directory("RESLT");
426
427 // Output number
428 doc_info.number()=0;
429
430 // Open a trace file
431 ofstream trace_file;
432 char filename[100];
433 sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
434 trace_file.open(filename);
435 trace_file << "VARIABLES=\"time\",\"u<SUB>FE</SUB>\","
436 << "\"u<SUB>exact</SUB>\",\"norm of error\",\"norm of solution\""
437 << std::endl;
438
439 // Choose simulation interval and timestep
440 double t_max=0.5;
441 double dt=0.01;
442
443 // Initialise timestep -- also sets the weights for all timesteppers
444 // in the problem.
445 problem.initialise_dt(dt);
446
447 // Set IC
448 problem.set_initial_condition();
449
450 //Output initial condition
451 problem.doc_solution(doc_info,trace_file);
452
453 //Increment counter for solutions
454 doc_info.number()++;
455
456 // Find number of steps
457 unsigned nstep = unsigned(t_max/dt);
458
459 // Timestepping loop
460 for (unsigned istep=0;istep<nstep;istep++)
461 {
462 cout << " Timestep " << istep << std::endl;
463
464 // Take timestep
465 problem.unsteady_newton_solve(dt);
466
467 //Output solution
468 problem.doc_solution(doc_info,trace_file);
469
470 //Increment counter for solutions
471 doc_info.number()++;
472 }
473
474 // Close trace file
475 trace_file.close();
476
477
478}; // end of main
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc the solution.
void actions_before_implicit_timestep()
Update the problem specs before next timestep: Set Dirchlet boundary conditions from exact solution.
void actions_after_newton_solve()
Update the problem specs after solve (empty)
UnsteadyHeatEquations< 2 >::UnsteadyHeatSourceFctPt Source_fct_pt
Pointer to source function.
void set_initial_condition()
Set initial condition (incl previous timesteps) according to specified function.
void actions_before_newton_solve()
Update the problem specs before solve (empty)
UnsteadyHeatProblem(UnsteadyHeatEquations< 2 >::UnsteadyHeatSourceFctPt source_fct_pt)
Constructor.
Node * Control_node_pt
Pointer to control node at which the solution is documented.
~UnsteadyHeatProblem()
Destructor (empty)
void actions_after_implicit_timestep()
Update the problem specs after solve (empty)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void get_exact_u(const double &time, const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
double Phi
Angle of bump.
void get_source(const double &time, const Vector< double > &x, double &source)
Source function to make it an exact solution.
int main()
///////////////////////////////////////////////////////////////////// ///////////////////////////////...