two_d_linear_wave.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 solving the 2D wave equation in a rectangle, with a
27// step wave moving through it at an angle.
28
29//Generic includes
30#include "generic.h"
31
32//Linear wave includes
33#include "linear_wave.h"
34
35// Mesh
36#include "meshes/rectangular_quadmesh.h"
37
38using namespace std;
39
40using namespace oomph;
41
42using namespace MathematicalConstants;
43
44//==start_of_tanh_solution============================================
45/// Namespace for exact solution for LinearWave equation
46/// with sharp step
47//====================================================================
49{
50
51 /// Parameter for steepness of step
52 double Alpha;
53
54 /// Orientation of step wave
55 double Phi;
56
57 /// Exact solution
58 double exact_u(const double& time, const Vector<double>& x)
59 {
60 double zeta=cos(Phi)*x[0]+sin(Phi)*x[1];
61 return tanh(1.0-Alpha*(zeta-time));
62 }
63
64 /// 1st time-deriv of exact solution
65 double exact_dudt(const double& time, const Vector<double>& x)
66 {
67 double zeta=cos(Phi)*x[0]+sin(Phi)*x[1];
68 return Alpha/(cosh(1.0-Alpha*(zeta-time))*
69 cosh(1.0-Alpha*(zeta-time)));
70 }
71
72 /// 2nd time-deriv of exact solution
73 double exact_d2udt2(const double& time, const Vector<double>& x)
74 {
75 double zeta=cos(Phi)*x[0]+sin(Phi)*x[1];
76 return -2.0*Alpha*Alpha*tanh(1.0-Alpha*(zeta-time))/
77 (cosh(1.0-Alpha*(zeta-time))*cosh(1.0-Alpha*(zeta-time)));
78 }
79
80
81 /// Exact solution as a vector
82 void get_exact_u(const double& time, const Vector<double>& x,
83 Vector<double>& u)
84 {
85 u[0]=exact_u(time,x);
86 u[1]=exact_dudt(time,x);
87 u[2]=exact_d2udt2(time,x);
88 }
89
90 /// Source function to make it an exact solution
91 void get_source(const double& time, const Vector<double>& x, double& source)
92 {
93 source=0.0;
94 }
95
96} // end of tanh solution
97
98
99
100/// /////////////////////////////////////////////////////////////////////
101/// /////////////////////////////////////////////////////////////////////
102/// /////////////////////////////////////////////////////////////////////
103
104//===start_of_problem_class===========================================
105/// LinearWave problem in rectanglular domain
106//====================================================================
107template<class ELEMENT, class TIMESTEPPER>
108class LinearWaveProblem : public Problem
109{
110
111public:
112
113 /// Constructor: pass number of elements in x and y directions,
114 /// bool indicating impulsive or "smooth" start,
115 /// and pointer to source function
116 LinearWaveProblem(const unsigned& nx, const unsigned& ny,
117 const bool& impulsive_start,
118 LinearWaveEquations<2>::LinearWaveSourceFctPt source_fct_pt);
119
120 /// Destructor (empty)
122
123 /// Update the problem specs after solve (empty)
125
126 /// Update the problem specs before solve (empty)
128
129 /// Update the problem specs after solve (empty)
131
132 /// Update the problem specs before next timestep:
133 /// Set time-dependent Dirchlet boundary from exact solution.
135 {
136 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
137 initial_value_fct(1);
138 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
139 initial_veloc_fct(1);
140 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
141 initial_accel_fct(1);
142
143 // Assign values for analytical value, veloc and accel:
144 initial_value_fct[0]=&TanhSolnForLinearWave::exact_u;
145 initial_veloc_fct[0]=&TanhSolnForLinearWave::exact_dudt;
146 initial_accel_fct[0]=&TanhSolnForLinearWave::exact_d2udt2;
147
148 // Loop over boundaries
149 unsigned num_bound=mesh_pt()->nboundary();
150 for (unsigned ibound=0;ibound<num_bound;ibound++)
151 {
152 // Loop over boundary nodes
153 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
154 for (unsigned inod=0;inod<num_nod;inod++)
155 {
156 // Set the boundary condition from the exact solution
157 Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
158
159 bool use_direct_assignment=false;
160 if (use_direct_assignment)
161 {
162 // Set nodal coordinates for evaluation of BC:
163 Vector<double> x(2);
164 x[0]=nod_pt->x(0);
165 x[1]=nod_pt->x(1);
166
167 // Set exact solution at current time
168 nod_pt->set_value(0,
169 TanhSolnForLinearWave::exact_u(time_pt()->time(),x));
170 }
171 else
172 {
173 // Get timestepper
174 TIMESTEPPER* timestepper_pt=dynamic_cast<TIMESTEPPER*>
175 (time_stepper_pt());
176
177 // Assign the history values
178 timestepper_pt->assign_initial_data_values(nod_pt,
179 initial_value_fct,
180 initial_veloc_fct,
181 initial_accel_fct);
182 }
183 }
184 }
185 } // end of actions before timestep
186
187 /// Set initial condition (incl history values)
189
190 /// Doc the solution
191 void doc_solution(DocInfo& doc_info);
192
193 /// Do unsteady run
194 void unsteady_run();
195
196private:
197
198 // Trace file
199 ofstream Trace_file;
200
201 // Impulsive start?
203
204}; // end of problem class
205
206
207
208//===start_of_constructor=================================================
209/// Constructor for LinearWave problem
210//========================================================================
211template<class ELEMENT, class TIMESTEPPER>
213 const unsigned& nx, const unsigned& ny, const bool& impulsive_start,
214 LinearWaveEquations<2>::LinearWaveSourceFctPt source_fct_pt) :
215 Impulsive_start(impulsive_start)
216{
217
218 //Allocate the timestepper -- this constructs the time object as well
219 add_time_stepper_pt(new TIMESTEPPER());
220
221 // Set up parameters for exact solution
222 //-------------------------------------
223
224 // Steepness of tanh profile
226
227 // Orientation of step wave
228 TanhSolnForLinearWave::Phi=MathematicalConstants::Pi/180.0*30.0;
229
230 // Set up mesh
231 //------------
232
233 // # of elements in x-direction
234 unsigned Nx=nx;
235
236 // # of elements in y-direction
237 unsigned Ny=ny;
238
239 // Domain length in x-direction
240 double Lx=1.0;
241
242 // Domain length in y-direction
243 double Ly=2.0;
244
245 // Build and assign mesh
246 Problem::mesh_pt()=new RectangularQuadMesh<ELEMENT>(
247 Nx,Ny,Lx,Ly,time_stepper_pt());
248
249 // Set the boundary conditions for this problem: All nodes are
250 // free by default -- just pin the ones that have Dirichlet conditions
251 // here.
252 unsigned num_bound = mesh_pt()->nboundary();
253 for(unsigned ibound=0;ibound<num_bound;ibound++)
254 {
255 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
256 for (unsigned inod=0;inod<num_nod;inod++)
257 {
258 mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
259 }
260 } //end of boundary conditions
261
262 // Complete build of elements
263 // --------------------------
264
265 // Loop over the elements to set up element-specific
266 // things that cannot be handled by constructor
267unsigned n_element = mesh_pt()->nelement();
268for(unsigned i=0;i<n_element;i++)
269 {
270 // Upcast from GeneralisedElement to the present element
271 ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
272
273 //Set the source function pointer
274 el_pt->source_fct_pt() = source_fct_pt;
275 }
276
277 // Setup equation numbering scheme
278 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
279
280} // end of constructor
281
282
283
284
285//===start_of_set_initial_condition=======================================
286/// Set initial condition.
287//========================================================================
288template<class ELEMENT, class TIMESTEPPER>
290{
291
292 // Get timestepper
293 TIMESTEPPER* timestepper_pt=dynamic_cast<TIMESTEPPER*>(time_stepper_pt());
294
295
296 // Impulsive start
297 //----------------
298 if (Impulsive_start)
299 {
300 // Loop over the nodes to set initial conditions everywhere
301 unsigned num_nod=mesh_pt()->nnode();
302 for (unsigned jnod=0;jnod<num_nod;jnod++)
303 {
304 // Pointer to node
305 Node* nod_pt=mesh_pt()->node_pt(jnod);
306
307 // Get nodal coordinates
308 Vector<double> x(2);
309 x[0]=nod_pt->x(0);
310 x[1]=nod_pt->x(1);
311
312 // Assign initial value from exact solution
313 nod_pt->set_value(0,TanhSolnForLinearWave::exact_u(time_pt()->time(),x));
314
315 // Set history values so that they are consistent with an impulsive
316 // start from this value
317 timestepper_pt->assign_initial_values_impulsive(nod_pt);
318 }
319 } // end impulsive start
320
321 // "Smooth" start from analytical time history
322 //--------------------------------------------
323 else
324 {
325
326 // Vector of function pointers to functions that specify the
327 // value, and the first and second time-derivatives of the
328 // function used as the initial condition
329 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
330 initial_value_fct(1);
331 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
332 initial_veloc_fct(1);
333 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
334 initial_accel_fct(1);
335
336 // Assign values for analytical value, veloc and accel:
337 initial_value_fct[0]=&TanhSolnForLinearWave::exact_u;
338 initial_veloc_fct[0]=&TanhSolnForLinearWave::exact_dudt;
339 initial_accel_fct[0]=&TanhSolnForLinearWave::exact_d2udt2;
340
341 // Assign Newmark history values so that Newmark approximations
342 // for velocity and accel are correct at initial time:
343
344 // Loop over the nodes to set initial conditions everywhere
345 unsigned num_nod=mesh_pt()->nnode();
346 for (unsigned jnod=0;jnod<num_nod;jnod++)
347 {
348 // Pointer to node
349 Node* nod_pt=mesh_pt()->node_pt(jnod);
350
351 // Assign the history values
352 timestepper_pt->assign_initial_data_values(nod_pt,
353 initial_value_fct,
354 initial_veloc_fct,
355 initial_accel_fct);
356 } // end of smooth start
357
358
359 // Paranoia: Check that the initial values were assigned correctly
360 double err_max=0.0;
361 for (unsigned jnod=0;jnod<num_nod;jnod++)
362 {
363 // Pointer to node
364 Node* nod_pt=mesh_pt()->node_pt(jnod);
365
366 // Get nodal coordinates
367 Vector<double> x(2);
368 x[0]=nod_pt->x(0);
369 x[1]=nod_pt->x(1);
370
371 // Get exact value and first and second time-derivatives
372 double u_exact=
373 TanhSolnForLinearWave::exact_u(time_pt()->time(),x);
374 double dudt_exact=
375 TanhSolnForLinearWave::exact_dudt(time_pt()->time(),x);
376 double d2udt2_exact=
377 TanhSolnForLinearWave::exact_d2udt2(time_pt()->time(),x);
378
379 // Get Newmark approximations for zero-th, first and second
380 // time-derivatives of the nodal values.
381 double u_fe=timestepper_pt->time_derivative(0,nod_pt,0);
382 double dudt_fe=timestepper_pt->time_derivative(1,nod_pt,0);
383 double d2udt2_fe=timestepper_pt->time_derivative(2,nod_pt,0);
384
385 // Error
386 double error=sqrt(pow(u_exact-u_fe,2)+
387 pow(dudt_exact-dudt_fe,2)+
388 pow(d2udt2_exact-d2udt2_fe,2));
389 if (error>err_max) err_max=error;
390 }
391 cout << "Max. error in assignment of initial condition "
392 << err_max << std::endl;
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, class TIMESTEPPER>
405{
406
407 ofstream some_file;
408 char filename[100];
409
410 // Number of plot points
411 unsigned npts;
412 npts=5;
413
414 cout << std::endl;
415 cout << "=================================================" << std::endl;
416 cout << "Docing solution for t=" << time_pt()->time() << std::endl;
417 cout << "=================================================" << std::endl;
418
419 // Output solution
420 //-----------------
421 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
422 doc_info.number());
423 some_file.open(filename);
424 mesh_pt()->output(some_file,npts);
425 some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
426 << time_pt()->time() << "\"";
427 some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
428 some_file << "1" << std::endl;
429 some_file << "2" << std::endl;
430 some_file << " 0 0" << std::endl;
431 some_file << time_pt()->time()*20.0 << " 0" << std::endl;
432 some_file.close();
433
434 // Output exact solution
435 //----------------------
436 sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),
437 doc_info.number());
438 oomph_info << " FILENAME: " << filename << std::endl;
439 some_file.open(filename);
440 mesh_pt()->output_fct(some_file,npts,time_pt()->time(),
442 some_file.close();
443
444 // Doc error
445 //----------
446 double error,norm;
447 sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),
448 doc_info.number());
449 some_file.open(filename);
450 mesh_pt()->compute_error(some_file,
452 time_pt()->time(),
453 error,norm);
454 some_file.close();
455 cout << "error: " << error << std::endl;
456 cout << "norm : " << norm << std::endl << std::endl;
457
458 // Write trace file
459 Trace_file << time_pt()->time() << " " << time_pt()->dt()
460 << " " << mesh_pt()->nelement() << " "
461 << error << " " << norm << std::endl;
462
463} // end of doc solution
464
465
466
467
468//===start_of_unsteady_run================================================
469/// Unsteady run.
470//========================================================================
471template<class ELEMENT, class TIMESTEPPER>
473{
474
475 // Setup labels for output
476 DocInfo doc_info;
477
478 // Output directory
479 if (Impulsive_start)
480 {
481 doc_info.set_directory("RESLT_impulsive");
482 }
483 else
484 {
485 doc_info.set_directory("RESLT_smooth");
486 }
487
488 // Open trace file
489 char filename[100];
490 sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
491 Trace_file.open(filename);
492
493 // Initialise time
494 double time0=0.0;
495 time_pt()->time()=time0;
496
497 // Set initial timestep
498 double dt=0.005;
499 time_pt()->initialise_dt(dt);
500
501 // Set IC
502 set_initial_condition();
503
504 //Output initial condition
505 doc_solution(doc_info);
506
507 //Increment counter for solutions
508 doc_info.number()++;
509
510 // Maximum time
511 double t_max=4.0;
512
513 // Number of steps
514 unsigned nstep=unsigned(t_max/dt);
515
516 // If validation run only do 2 timesteps
517 if (CommandLineArgs::Argc>1)
518 {
519 nstep=2;
520 cout << "Validation run -- only doing two timesteps." << std::endl;
521 }
522
523 // Timestepping loop
524 for (unsigned istep=0;istep<nstep;istep++)
525 {
526 //Take fixed timestep without spatial adaptivity
527 unsteady_newton_solve(dt);
528
529 //Output solution
530 doc_solution(doc_info);
531
532 //Increment counter for solutions
533 doc_info.number()++;
534 }
535
536 // Close trace file
537 Trace_file.close();
538
539} // end of unsteady run
540
541
542
543
544/// /////////////////////////////////////////////////////////////////////
545/// /////////////////////////////////////////////////////////////////////
546/// /////////////////////////////////////////////////////////////////////
547
548
549
550//===start_of_main========================================================
551/// Demonstrate how to solve LinearWave problem.
552//========================================================================
553int main(int argc, char* argv[])
554{
555
556 // Store command line arguments: If a command line argument is specied
557 // we regard this as validation run.
558 CommandLineArgs::setup(argc,argv);
559
560 // Number of elements in x direction
561 unsigned n_x=10;
562
563 // Number of elements in y direction
564 unsigned n_y=20;
565
566 // Impulsive start?
567 bool impulsive_start;
568
569 // Run with impulsive start
570 // ------------------------
571 {
572 impulsive_start=true;
573
574 // Build problem
576 problem(n_x,n_y,impulsive_start,&TanhSolnForLinearWave::get_source);
577
578 // Run it
579 problem.unsteady_run();
580 }
581
582 // Run with "smooth" start
583 // -----------------------
584 {
585 impulsive_start=false;
586
587 // Build problem
589 problem(n_x,n_y,impulsive_start,&TanhSolnForLinearWave::get_source);
590
591 // Run it
592 problem.unsteady_run();
593 }
594
595
596}; // end of main
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void actions_before_implicit_timestep()
Update the problem specs before next timestep: Set time-dependent Dirchlet boundary from exact soluti...
LinearWaveProblem(const unsigned &nx, const unsigned &ny, const bool &impulsive_start, LinearWaveEquations< 2 >::LinearWaveSourceFctPt source_fct_pt)
Constructor: pass number of elements in x and y directions, bool indicating impulsive or "smooth" sta...
void actions_after_implicit_timestep()
Update the problem specs after solve (empty)
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void actions_before_newton_solve()
Update the problem specs before solve (empty)
void doc_solution(DocInfo &doc_info)
Doc the solution.
void unsteady_run()
Do unsteady run.
~LinearWaveProblem()
Destructor (empty)
void set_initial_condition()
Set initial condition (incl history values)
Namespace for exact solution for LinearWave equation with sharp step.
void get_source(const double &time, const Vector< double > &x, double &source)
Source function to make it an exact solution.
double Alpha
Parameter for steepness of step.
double Phi
Orientation of step wave.
double exact_d2udt2(const double &time, const Vector< double > &x)
2nd time-deriv of exact solution
void get_exact_u(const double &time, const Vector< double > &x, Vector< double > &u)
Exact solution as a vector.
double exact_dudt(const double &time, const Vector< double > &x)
1st time-deriv of exact solution
double exact_u(const double &time, const Vector< double > &x)
Exact solution.
int main(int argc, char *argv[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...