32 #include "unsteady_heat.h"
35 #include "meshes/rectangular_quadmesh.h"
39 using namespace oomph;
41 using namespace MathematicalConstants;
61 void get_exact_u(
const double& time,
const Vector<double>& x,
64 double zeta=cos(
Phi)*x[0]+sin(
Phi)*x[1];
65 u[0]=exp(-
K*time)*sin(zeta*sqrt(
K));
69 void get_exact_u(
const double& time,
const Vector<double>& x,
double& u)
71 double zeta=cos(
Phi)*x[0]+sin(
Phi)*x[1];
72 u=exp(-
K*time)*sin(zeta*sqrt(
K));
76 void get_source(
const double& time,
const Vector<double>& x,
double& source)
90 template<
class ELEMENT>
124 void dump_it(ofstream& dump_file);
127 void restart(ifstream& restart_file);
132 UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt Source_fct_pt;
135 Node* Control_node_pt;
143 template<
class ELEMENT>
145 UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt source_fct_pt) :
146 Source_fct_pt(source_fct_pt)
152 add_time_stepper_pt(
new BDF<2>);
172 mesh_pt() =
new RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,time_stepper_pt());
177 unsigned n_el=mesh_pt()->nelement();
180 unsigned control_el=unsigned(n_el/2);
185 cout <<
"Recording trace of the solution at: "
194 unsigned n_bound = mesh_pt()->nboundary();
195 for(
unsigned b=0;b<n_bound;b++)
197 unsigned n_node = mesh_pt()->nboundary_node(b);
198 for (
unsigned n=0;n<n_node;n++)
200 mesh_pt()->boundary_node_pt(b,n)->pin(0);
209 unsigned n_element = mesh_pt()->nelement();
213 for(
unsigned i=0;i<n_element;i++)
216 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(mesh_pt()->element_pt(i));
223 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
233 template<
class ELEMENT>
237 double time=time_pt()->time();
240 unsigned num_bound = mesh_pt()->nboundary();
241 for(
unsigned ibound=0;ibound<num_bound;ibound++)
244 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
245 for (
unsigned inod=0;inod<num_nod;inod++)
247 Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
255 nod_pt->set_value(0,u);
266 template<
class ELEMENT>
271 ifstream* restart_file_pt=0;
277 if (CommandLineArgs::Argc==1)
279 cout <<
"No restart -- setting IC from exact solution" << std::endl;
281 else if (CommandLineArgs::Argc==2)
284 restart_file_pt=
new ifstream(CommandLineArgs::Argv[1],ios_base::in);
285 if (restart_file_pt!=0)
287 cout <<
"Have opened " << CommandLineArgs::Argv[1] <<
288 " for restart. " << std::endl;
292 std::ostringstream error_stream;
294 <<
"ERROR while trying to open " << CommandLineArgs::Argv[1] <<
295 " for restart." << std::endl;
299 OOMPH_CURRENT_FUNCTION,
300 OOMPH_EXCEPTION_LOCATION);
306 std::ostringstream error_stream;
307 error_stream <<
"Can only specify one input file\n"
308 <<
"You specified the following command line arguments:\n";
310 CommandLineArgs::output();
314 OOMPH_CURRENT_FUNCTION,
315 OOMPH_EXCEPTION_LOCATION);
321 if (restart_file_pt!=0)
325 restart(*restart_file_pt);
341 double backed_up_time=time_pt()->time();
348 Vector<double> soln(1);
352 unsigned num_nod = mesh_pt()->nnode();
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--)
359 prev_time[itime]=time_pt()->time(
unsigned(itime));
363 for (
int itime=nprev_steps;itime>=0;itime--)
366 double time=prev_time[itime];
367 cout <<
"setting IC at time =" << time << std::endl;
370 for (
unsigned n=0;n<num_nod;n++)
373 x[0]=mesh_pt()->node_pt(n)->x(0);
374 x[1]=mesh_pt()->node_pt(n)->x(1);
380 mesh_pt()->node_pt(n)->set_value(itime,0,soln[0]);
383 for (
unsigned i=0;i<2;i++)
385 mesh_pt()->node_pt(n)->x(itime,i)=x[i];
391 time_pt()->time()=backed_up_time;
403 template<
class ELEMENT>
416 cout <<
"=================================================" << std::endl;
417 cout <<
"Docing solution for t=" << time_pt()->time() << std::endl;
418 cout <<
"=================================================" << std::endl;
423 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
425 some_file.open(filename);
426 mesh_pt()->output(some_file,npts);
429 some_file <<
"TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
430 << time_pt()->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;
443 sprintf(filename,
"%s/exact_soln%i.dat",doc_info.directory().c_str(),
445 some_file.open(filename);
446 mesh_pt()->output_fct(some_file,npts,time_pt()->time(),
453 sprintf(filename,
"%s/error%i.dat",doc_info.directory().c_str(),
455 some_file.open(filename);
456 mesh_pt()->compute_error(some_file,
464 cout <<
"error: " << error << std::endl;
465 cout <<
"norm : " << norm << std::endl << std::endl;
468 Vector<double> x_ctrl(2);
469 x_ctrl[0]=Control_node_pt->x(0);
470 x_ctrl[1]=Control_node_pt->x(1);
473 trace_file << time_pt()->time() <<
" "
474 << Control_node_pt->value(0) <<
" "
477 << norm << std::endl;
481 sprintf(filename,
"%s/restart%i.dat",doc_info.directory().c_str(),
483 some_file.open(filename);
495 template<
class ELEMENT>
500 Problem::dump(dump_file);
509 template<
class ELEMENT>
514 Problem::read(restart_file);
532 int main(
int argc,
char* argv[])
536 CommandLineArgs::setup(argc,argv);
546 doc_info.set_directory(
"RESLT");
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\""
578 double dt=problem.time_pt()->dt();
581 unsigned nstep = unsigned(t_max/dt);
584 for (
unsigned istep=0;istep<nstep;istep++)
587 problem.unsteady_newton_solve(dt);
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
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.
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[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...