two_d_unsteady_heat_restarted.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 with optional restart from disk
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
123 /// Dump problem to disk to allow for restart.
124 void dump_it(ofstream& dump_file);
125
126 /// Read problem for restart from specified restart file.
127 void restart(ifstream& restart_file);
128
129private:
130
131 /// Pointer to source function
132 UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt Source_fct_pt;
133
134 /// Pointer to control node at which the solution is documented
135 Node* Control_node_pt;
136
137}; // end of problem class
138
139
140//========start_of_constructor============================================
141/// Constructor for UnsteadyHeat problem in square domain
142//========================================================================
143template<class ELEMENT>
145 UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt source_fct_pt) :
146 Source_fct_pt(source_fct_pt)
147{
148
149 // Allocate the timestepper -- this constructs the Problem's
150 // time object with a sufficient amount of storage to store the
151 // previous timsteps.
152 add_time_stepper_pt(new BDF<2>);
153
154 // Setup parameters for exact solution
155 // -----------------------------------
156
157 // Decay parameter
159
160 // Setup mesh
161 //-----------
162
163 // Number of elements in x and y directions
164 unsigned nx=5;
165 unsigned ny=5;
166
167 // Lengths in x and y directions
168 double lx=1.0;
169 double ly=1.0;
170
171 // Build mesh
172 mesh_pt() = new RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,time_stepper_pt());
173
174 // Choose a control node at which the solution is documented
175 //----------------------------------------------------------
176 // Total number of elements
177 unsigned n_el=mesh_pt()->nelement();
178
179 // Choose an element in the middle
180 unsigned control_el=unsigned(n_el/2);
181
182 // Choose its first node as the control node
183 Control_node_pt=mesh_pt()->finite_element_pt(control_el)->node_pt(0);
184
185 cout << "Recording trace of the solution at: "
186 << Control_node_pt->x(0) << " "
187 << Control_node_pt->x(1) << std::endl;
188
189
190 // Set the boundary conditions for this problem:
191 // ---------------------------------------------
192 // All nodes are free by default -- just pin the ones that have
193 // Dirichlet conditions here.
194 unsigned n_bound = mesh_pt()->nboundary();
195 for(unsigned b=0;b<n_bound;b++)
196 {
197 unsigned n_node = mesh_pt()->nboundary_node(b);
198 for (unsigned n=0;n<n_node;n++)
199 {
200 mesh_pt()->boundary_node_pt(b,n)->pin(0);
201 }
202 } // end of set boundary conditions
203
204
205 // Complete the build of all elements so they are fully functional
206 //----------------------------------------------------------------
207
208 // Find number of elements in mesh
209 unsigned n_element = mesh_pt()->nelement();
210
211 // Loop over the elements to set up element-specific
212 // things that cannot be handled by constructor
213 for(unsigned i=0;i<n_element;i++)
214 {
215 // Upcast from FiniteElement 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
222 // Do equation numbering
223 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
224
225} // end of constructor
226
227
228
229//=========start of actions_before_implicit_timestep======================
230/// Actions before timestep: update the domain, then reset the
231/// boundary conditions for the current time.
232//========================================================================
233template<class ELEMENT>
235{
236 // Get current time
237 double time=time_pt()->time();
238
239 //Loop over the boundaries
240 unsigned num_bound = mesh_pt()->nboundary();
241 for(unsigned ibound=0;ibound<num_bound;ibound++)
242 {
243 // Loop over the nodes on boundary
244 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
245 for (unsigned inod=0;inod<num_nod;inod++)
246 {
247 Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
248 double u;
249 Vector<double> x(2);
250 x[0]=nod_pt->x(0);
251 x[1]=nod_pt->x(1);
252 // Get current values of the boundary conditions from the
253 // exact solution
255 nod_pt->set_value(0,u);
256 }
257 }
258} // end of actions_before_implicit_timestep
259
260
261
262//======================start_of_set_initial_condition====================
263/// Set initial condition: Assign previous and current values
264/// from exact solution or from restart file.
265//========================================================================
266template<class ELEMENT>
268{
269
270 // Pointer to restart file
271 ifstream* restart_file_pt=0;
272
273 // Restart?
274 //---------
275 // Restart file specified via command line [all programs have at least
276 // a single command line argument: their name. Ignore this here.]
277 if (CommandLineArgs::Argc==1)
278 {
279 cout << "No restart -- setting IC from exact solution" << std::endl;
280 }
281 else if (CommandLineArgs::Argc==2)
282 {
283 // Open restart file
284 restart_file_pt= new ifstream(CommandLineArgs::Argv[1],ios_base::in);
285 if (restart_file_pt!=0)
286 {
287 cout << "Have opened " << CommandLineArgs::Argv[1] <<
288 " for restart. " << std::endl;
289 }
290 else
291 {
292 std::ostringstream error_stream;
293 error_stream
294 << "ERROR while trying to open " << CommandLineArgs::Argv[1] <<
295 " for restart." << std::endl;
296
297 throw OomphLibError(
298 error_stream.str(),
299 OOMPH_CURRENT_FUNCTION,
300 OOMPH_EXCEPTION_LOCATION);
301 }
302 }
303 // More than one command line argument?
304 else
305 {
306 std::ostringstream error_stream;
307 error_stream << "Can only specify one input file\n"
308 << "You specified the following command line arguments:\n";
309 //Fix this
310 CommandLineArgs::output();
311
312 throw OomphLibError(
313 error_stream.str(),
314 OOMPH_CURRENT_FUNCTION,
315 OOMPH_EXCEPTION_LOCATION);
316 }
317
318
319 // Read restart data:
320 //-------------------
321 if (restart_file_pt!=0)
322 {
323
324 // Read the problem data from the restart file
325 restart(*restart_file_pt);
326
327 }
328 // Assign initial condition from exact solution
329 //---------------------------------------------
330 else
331 {
332
333 // Choose timestep
334 double dt=0.01;
335
336 // Initialise timestep -- also sets the weights for all timesteppers
337 // in the problem.
338 initialise_dt(dt);
339
340 // Backup time in global Time object
341 double backed_up_time=time_pt()->time();
342
343 // Past history fo needs to be established for t=time0-deltat, ...
344 // Then provide current values (at t=time0) which will also form
345 // the initial guess for the first solve at t=time0+deltat
346
347 // Vector of exact solution value
348 Vector<double> soln(1);
349 Vector<double> x(2);
350
351 //Find number of nodes in mesh
352 unsigned num_nod = mesh_pt()->nnode();
353
354 // Set continuous times at previous timesteps
355 int nprev_steps=time_stepper_pt()->nprev_values();
356 Vector<double> prev_time(nprev_steps+1);
357 for (int itime=nprev_steps;itime>=0;itime--)
358 {
359 prev_time[itime]=time_pt()->time(unsigned(itime));
360 }
361
362 // Loop over current & previous timesteps
363 for (int itime=nprev_steps;itime>=0;itime--)
364 {
365 // Continuous time
366 double time=prev_time[itime];
367 cout << "setting IC at time =" << time << std::endl;
368
369 // Loop over the nodes to set initial guess everywhere
370 for (unsigned n=0;n<num_nod;n++)
371 {
372 // Get nodal coordinates
373 x[0]=mesh_pt()->node_pt(n)->x(0);
374 x[1]=mesh_pt()->node_pt(n)->x(1);
375
376 // Get intial solution
378
379 // Assign solution
380 mesh_pt()->node_pt(n)->set_value(itime,0,soln[0]);
381
382 // Loop over coordinate directions: Previous position = present position
383 for (unsigned i=0;i<2;i++)
384 {
385 mesh_pt()->node_pt(n)->x(itime,i)=x[i];
386 }
387 }
388 }
389
390 // Reset backed up time for global timestepper
391 time_pt()->time()=backed_up_time;
392
393 }
394
395
396} // end of set_initial_condition
397
398
399
400//=======start_of_doc_solution============================================
401/// Doc the solution
402//========================================================================
403template<class ELEMENT>
405doc_solution(DocInfo& doc_info,ofstream& trace_file)
406{
407 ofstream some_file;
408 char filename[100];
409
410 // Number of plot points
411 unsigned npts;
412 npts=5;
413
414
415 cout << std::endl;
416 cout << "=================================================" << std::endl;
417 cout << "Docing solution for t=" << time_pt()->time() << std::endl;
418 cout << "=================================================" << std::endl;
419
420
421 // Output solution
422 //-----------------
423 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
424 doc_info.number());
425 some_file.open(filename);
426 mesh_pt()->output(some_file,npts);
427
428 // Write file as a tecplot text object
429 some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
430 << time_pt()->time() << "\"";
431 // ...and draw a horizontal line whose length is proportional
432 // to the elapsed time
433 some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
434 some_file << "1" << std::endl;
435 some_file << "2" << std::endl;
436 some_file << " 0 0" << std::endl;
437 some_file << time_pt()->time()*20.0 << " 0" << std::endl;
438 some_file.close();
439
440
441 // Output exact solution
442 //----------------------
443 sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
444 doc_info.number());
445 some_file.open(filename);
446 mesh_pt()->output_fct(some_file,npts,time_pt()->time(),
448 some_file.close();
449
450 // Doc error
451 //----------
452 double error,norm;
453 sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
454 doc_info.number());
455 some_file.open(filename);
456 mesh_pt()->compute_error(some_file,
458 time_pt()->time(),
459 error,norm);
460 some_file.close();
461
462 // Doc solution and error
463 //-----------------------
464 cout << "error: " << error << std::endl;
465 cout << "norm : " << norm << std::endl << std::endl;
466
467 // Get exact solution at control node
468 Vector<double> x_ctrl(2);
469 x_ctrl[0]=Control_node_pt->x(0);
470 x_ctrl[1]=Control_node_pt->x(1);
471 double u_exact;
472 ExactSolnForUnsteadyHeat::get_exact_u(time_pt()->time(),x_ctrl,u_exact);
473 trace_file << time_pt()->time() << " "
474 << Control_node_pt->value(0) << " "
475 << u_exact << " "
476 << error << " "
477 << norm << std::endl;
478
479
480 // Write restart file
481 sprintf(filename,"%s/restart%i.dat",doc_info.directory().c_str(),
482 doc_info.number());
483 some_file.open(filename);
484 dump_it(some_file);
485 some_file.close();
486
487} // end of doc_solution
488
489
490
491
492//=====start_of_dump_it===================================================
493/// Dump the solution to disk to allow for restart
494//========================================================================
495template<class ELEMENT>
497{
498
499 // Call generic dump()
500 Problem::dump(dump_file);
501
502} // end of dump_it
503
504
505
506//=================start_of_restart=======================================
507/// Read solution from disk for restart
508//========================================================================
509template<class ELEMENT>
510void UnsteadyHeatProblem<ELEMENT>::restart(ifstream& restart_file)
511{
512
513 // Read the generic problem data from restart file
514 Problem::read(restart_file);
515
516} // end of restart
517
518
519
520
521/// /////////////////////////////////////////////////////////////////////
522/// /////////////////////////////////////////////////////////////////////
523/// /////////////////////////////////////////////////////////////////////
524
525
526
527//=======start_of_main====================================================
528/// Driver code for unsteady heat equation with option for
529/// restart from disk: Only a single command line argument is allowed.
530/// If specified it is interpreted as the name of the restart file.
531//========================================================================
532int main(int argc, char* argv[])
533{
534
535 // Store command line arguments
536 CommandLineArgs::setup(argc,argv);
537
538 // Build problem
541
542 // Setup labels for output
543 DocInfo doc_info;
544
545 // Output directory
546 doc_info.set_directory("RESLT");
547
548 // Output number
549 doc_info.number()=0;
550
551 // Open a trace file
552 ofstream trace_file;
553 char filename[100];
554 sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
555 trace_file.open(filename);
556 trace_file << "VARIABLES=\"time\",\"u<SUB>FE</SUB>\","
557 << "\"u<SUB>exact</SUB>\",\"norm of error\",\"norm of solution\""
558 << std::endl;
559
560 // Choose simulation interval
561 double t_max=0.5;
562
563 // Set IC
564 problem.set_initial_condition();
565
566 //Output initial condition
567 problem.doc_solution(doc_info,trace_file);
568
569 //Increment counter for solutions
570 doc_info.number()++;
571
572
573 // Get (fixed) timestep dt from Time object. dt was set in
574 // set_initial_conditions(), either by assigning a
575 // value or re-cycling the one from the restart file.
576 // This is necessary because the history values obtained from the
577 // restart file are only correct if the same dts are used.
578 double dt=problem.time_pt()->dt();
579
580 // Find number of steps
581 unsigned nstep = unsigned(t_max/dt);
582
583 // Timestepping loop
584 for (unsigned istep=0;istep<nstep;istep++)
585 {
586 // Take timestep
587 problem.unsteady_newton_solve(dt);
588
589 //Output solution
590 problem.doc_solution(doc_info,trace_file);
591
592 //Increment counter for solutions
593 doc_info.number()++;
594 }
595
596 // Close trace file
597 trace_file.close();
598
599
600}; // 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 dump_it(ofstream &dump_file)
Dump problem to disk to allow for restart.
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 restart(ifstream &restart_file)
Read problem for restart from specified restart file.
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(int argc, char *argv[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...