33 #include "linear_wave.h"
36 #include "meshes/rectangular_quadmesh.h"
40 using namespace oomph;
42 using namespace MathematicalConstants;
58 double exact_u(
const double& time,
const Vector<double>& x)
60 double zeta=cos(
Phi)*x[0]+sin(
Phi)*x[1];
61 return tanh(1.0-
Alpha*(zeta-time));
65 double exact_dudt(
const double& time,
const Vector<double>& x)
67 double zeta=cos(
Phi)*x[0]+sin(
Phi)*x[1];
69 cosh(1.0-
Alpha*(zeta-time)));
75 double zeta=cos(
Phi)*x[0]+sin(
Phi)*x[1];
77 (cosh(1.0-
Alpha*(zeta-time))*cosh(1.0-
Alpha*(zeta-time)));
91 void get_source(
const double& time,
const Vector<double>& x,
double& source)
107 template<
class ELEMENT,
class TIMESTEPPER>
117 const bool& impulsive_start,
118 LinearWaveEquations<2>::LinearWaveSourceFctPt source_fct_pt);
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);
149 unsigned num_bound=mesh_pt()->nboundary();
150 for (
unsigned ibound=0;ibound<num_bound;ibound++)
153 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
154 for (
unsigned inod=0;inod<num_nod;inod++)
157 Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
159 bool use_direct_assignment=
false;
160 if (use_direct_assignment)
174 TIMESTEPPER* timestepper_pt=
dynamic_cast<TIMESTEPPER*
>
178 timestepper_pt->assign_initial_data_values(nod_pt,
188 void set_initial_condition();
191 void doc_solution(DocInfo& doc_info);
211 template<
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)
219 add_time_stepper_pt(
new TIMESTEPPER());
246 Problem::mesh_pt()=
new RectangularQuadMesh<ELEMENT>(
247 Nx,Ny,Lx,Ly,time_stepper_pt());
252 unsigned num_bound = mesh_pt()->nboundary();
253 for(
unsigned ibound=0;ibound<num_bound;ibound++)
255 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
256 for (
unsigned inod=0;inod<num_nod;inod++)
258 mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
267 unsigned n_element = mesh_pt()->nelement();
268 for(
unsigned i=0;i<n_element;i++)
271 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(mesh_pt()->element_pt(i));
274 el_pt->source_fct_pt() = source_fct_pt;
278 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
288 template<
class ELEMENT,
class TIMESTEPPER>
293 TIMESTEPPER* timestepper_pt=
dynamic_cast<TIMESTEPPER*
>(time_stepper_pt());
301 unsigned num_nod=mesh_pt()->nnode();
302 for (
unsigned jnod=0;jnod<num_nod;jnod++)
305 Node* nod_pt=mesh_pt()->node_pt(jnod);
317 timestepper_pt->assign_initial_values_impulsive(nod_pt);
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);
345 unsigned num_nod=mesh_pt()->nnode();
346 for (
unsigned jnod=0;jnod<num_nod;jnod++)
349 Node* nod_pt=mesh_pt()->node_pt(jnod);
352 timestepper_pt->assign_initial_data_values(nod_pt,
361 for (
unsigned jnod=0;jnod<num_nod;jnod++)
364 Node* nod_pt=mesh_pt()->node_pt(jnod);
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);
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;
391 cout <<
"Max. error in assignment of initial condition "
392 << err_max << std::endl;
403 template<
class ELEMENT,
class TIMESTEPPER>
415 cout <<
"=================================================" << std::endl;
416 cout <<
"Docing solution for t=" << time_pt()->time() << std::endl;
417 cout <<
"=================================================" << std::endl;
421 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
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;
436 sprintf(filename,
"%s/exact_soln%i.dat",doc_info.directory().c_str(),
438 oomph_info <<
" FILENAME: " << filename << std::endl;
439 some_file.open(filename);
440 mesh_pt()->output_fct(some_file,npts,time_pt()->time(),
447 sprintf(filename,
"%s/error%i.dat",doc_info.directory().c_str(),
449 some_file.open(filename);
450 mesh_pt()->compute_error(some_file,
455 cout <<
"error: " << error << std::endl;
456 cout <<
"norm : " << norm << std::endl << std::endl;
459 Trace_file << time_pt()->time() <<
" " << time_pt()->dt()
460 <<
" " << mesh_pt()->nelement() <<
" "
461 << error <<
" " << norm << std::endl;
471 template<
class ELEMENT,
class TIMESTEPPER>
481 doc_info.set_directory(
"RESLT_impulsive");
485 doc_info.set_directory(
"RESLT_smooth");
490 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
491 Trace_file.open(filename);
495 time_pt()->time()=time0;
499 time_pt()->initialise_dt(dt);
502 set_initial_condition();
505 doc_solution(doc_info);
514 unsigned nstep=unsigned(t_max/dt);
517 if (CommandLineArgs::Argc>1)
520 cout <<
"Validation run -- only doing two timesteps." << std::endl;
524 for (
unsigned istep=0;istep<nstep;istep++)
527 unsteady_newton_solve(dt);
530 doc_solution(doc_info);
553 int main(
int argc,
char* argv[])
558 CommandLineArgs::setup(argc,argv);
567 bool impulsive_start;
572 impulsive_start=
true;
585 impulsive_start=
false;
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
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[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...