32 #include "unsteady_heat.h"
35 #include "meshes/quarter_circle_sector_mesh.h"
39 using namespace oomph;
41 using namespace MathematicalConstants;
57 const double& a_hat,
const double& b_hat,
58 const double& period, Time* time_pt) :
59 GeomObject(1,2), A(a), B(b), A_hat(a_hat), B_hat(b_hat),
60 T(period), Time_pt(time_pt) {}
67 void position(
const Vector<double>& xi, Vector<double>& r)
const
70 double time=Time_pt->time();
73 r[0] = (A+A_hat*sin(2.0*MathematicalConstants::Pi*time/T))*cos(xi[0]);
74 r[1] = (B+B_hat*sin(2.0*MathematicalConstants::Pi*time/T))*sin(xi[0]);
83 void position(
const unsigned& t,
const Vector<double>& xi,
84 Vector<double>& r)
const
87 double time=Time_pt->time(t);
90 r[0] = (A+A_hat*sin(2.0*MathematicalConstants::Pi*time/T))*cos(xi[0]);
91 r[1] = (B+B_hat*sin(2.0*MathematicalConstants::Pi*time/T))*sin(xi[0]);
145 return Beta*tanh(
Gamma*cos(0.2E1*MathematicalConstants::Pi*time));
155 MathematicalConstants::Pi*time)))-Y));
159 void get_exact_u(
const double& time,
const Vector<double>& x,
double& u)
164 MathematicalConstants::Pi*time)))-Y));
168 void get_source(
const double& time,
const Vector<double>& x,
double& source)
173 MathematicalConstants::Pi*time)))-Y))*(1.0-pow(tanh(0.1E1-
Alpha*(
TanPhi*(X-
174 Beta*tanh(
Gamma*cos(0.2E1*MathematicalConstants::Pi*time)))-Y)),2.0))*
Alpha*
176 MathematicalConstants::Pi*time)))-Y))*(1.0-pow(tanh(0.1E1-
Alpha*(
TanPhi*(X-
177 Beta*tanh(
Gamma*cos(0.2E1*MathematicalConstants::Pi*time)))-Y)),2.0))*
Alpha*
179 MathematicalConstants::Pi*time)))-Y)),2.0))*
Alpha*
TanPhi*
Beta*(1.0-pow(tanh(
180 Gamma*cos(0.2E1*MathematicalConstants::Pi*time)),2.0))*
Gamma*sin(0.2E1*
181 MathematicalConstants::Pi*time)*MathematicalConstants::Pi;
186 const Vector<double>& x,
198 MathematicalConstants::Pi*time)))-Y)),2.0))*
Alpha*
TanPhi*NX+(1.0-pow(tanh(
200 time)))-Y)),2.0))*
Alpha*NY;
214 template<
class ELEMENT>
222 UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt source_fct_pt);
238 void actions_before_implicit_timestep();
241 void actions_before_adapt();
244 void actions_after_adapt();
247 double global_temporal_error_norm();
251 void set_initial_condition();
256 void create_flux_elements(
const unsigned &b, Mesh*
const &bulk_mesh_pt,
257 Mesh*
const &surface_mesh_pt);
260 void delete_flux_elements(Mesh*
const &surface_mesh_pt);
266 void dump_it(ofstream& dump_file);
269 void restart(ifstream& restart_file);
281 void write_trace_file_header();
323 template<
class ELEMENT>
325 UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt source_fct_pt)
326 : Source_fct_pt(source_fct_pt)
344 add_time_stepper_pt(
new BDF<2>(
true));
386 double xi_hi=MathematicalConstants::Pi/2.0;
391 double fract_mid=0.5;
392 Bulk_mesh_pt =
new RefineableQuarterCircleSectorMesh<ELEMENT>(
393 Boundary_pt,xi_lo,fract_mid,xi_hi,time_stepper_pt());
411 Z2ErrorEstimator* error_estimator_pt=
new Z2ErrorEstimator;
412 Bulk_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
418 for(
unsigned b=0;b<n_bound;b++)
425 for (
unsigned n=0;n<n_node;n++)
434 FiniteElement* el0_pt=
Bulk_mesh_pt->finite_element_pt(0);
435 unsigned nnod=el0_pt->nnode();
447 for(
unsigned i=0;i<n_element;i++)
450 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(
Bulk_mesh_pt->element_pt(i));
458 for(
unsigned e=0;e<n_element;e++)
461 UnsteadyHeatFluxElement<ELEMENT> *el_pt =
462 dynamic_cast<UnsteadyHeatFluxElement<ELEMENT>*
>(
466 el_pt->flux_fct_pt() =
471 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
479 template<
class ELEMENT>
491 template<
class ELEMENT>
495 Bulk_mesh_pt->node_update();
498 double time=time_pt()->time();
501 unsigned num_bound = Bulk_mesh_pt->nboundary();
502 for(
unsigned b=0;b<num_bound;b++)
508 unsigned num_nod=Bulk_mesh_pt->nboundary_node(b);
509 for (
unsigned j=0;j<num_nod;j++)
511 Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
517 nod_pt->set_value(0,u);
527 template<
class ELEMENT>
532 delete_flux_elements(Surface_mesh_pt);
535 rebuild_global_mesh();
543 template<
class ELEMENT>
548 create_flux_elements(0,Bulk_mesh_pt,Surface_mesh_pt);
551 rebuild_global_mesh();
554 unsigned n_element=Surface_mesh_pt->nelement();
555 for(
unsigned e=0;e<n_element;e++)
558 UnsteadyHeatFluxElement<ELEMENT> *el_pt =
559 dynamic_cast<UnsteadyHeatFluxElement<ELEMENT>*
>(
560 Surface_mesh_pt->element_pt(e));
563 el_pt->flux_fct_pt() =
573 template<
class ELEMENT>
578 ifstream* restart_file_pt=0;
584 if (CommandLineArgs::Argc==1)
586 cout <<
"No restart -- setting IC from exact solution" << std::endl;
588 else if (CommandLineArgs::Argc==2)
591 restart_file_pt=
new ifstream(CommandLineArgs::Argv[1],ios_base::in);
592 if (restart_file_pt!=0)
594 cout <<
"Have opened " << CommandLineArgs::Argv[1] <<
595 " for restart. " << std::endl;
599 std::ostringstream error_stream;
601 <<
"ERROR while trying to open " << CommandLineArgs::Argv[1]
602 <<
" for restart." << std::endl;
606 OOMPH_CURRENT_FUNCTION,
607 OOMPH_EXCEPTION_LOCATION);
613 std::ostringstream error_stream;
614 error_stream <<
"Can only specify one input file\n"
615 <<
"You specified the following command line arguments:\n";
619 OOMPH_CURRENT_FUNCTION,
620 OOMPH_EXCEPTION_LOCATION);
624 CommandLineArgs::output();
630 if (restart_file_pt!=0)
634 restart(*restart_file_pt);
651 double backed_up_time=time_pt()->time();
658 Vector<double> soln(1);
662 unsigned num_nod = Bulk_mesh_pt->nnode();
665 int nprev_steps=time_stepper_pt()->nprev_values();
666 Vector<double> prev_time(nprev_steps+1);
667 for (
int itime=nprev_steps;itime>=0;itime--)
669 prev_time[itime]=time_pt()->time(
unsigned(itime));
674 for (
int itime=nprev_steps;itime>=0;itime--)
676 double time=prev_time[itime];
680 time_pt()->time()=time;
682 cout <<
"setting IC at time =" << time << std::endl;
687 Bulk_mesh_pt->node_update();
690 for (
unsigned jnod=0;jnod<num_nod;jnod++)
693 x[0]=Bulk_mesh_pt->node_pt(jnod)->x(0);
694 x[1]=Bulk_mesh_pt->node_pt(jnod)->x(1);
700 Bulk_mesh_pt->node_pt(jnod)->set_value(itime,0,soln[0]);
703 for (
unsigned i=0;i<2;i++)
705 Bulk_mesh_pt->node_pt(jnod)->x(itime,i)=x[i];
711 time_pt()->time()=backed_up_time;
722 template<
class ELEMENT>
727 sprintf(filename,
"%s/trace.dat",Doc_info.directory().c_str());
728 Trace_file.open(filename);
730 Trace_file <<
"VARIABLES=\"time t\","
731 <<
"\"u<SUB>FE</SUB>\","
732 <<
"\"u<SUB>exact</SUB>\","
735 <<
"\"global temporal error norm\","
736 <<
"\"X<SUB>step</SUB>\","
737 <<
"\"N<SUB>element</SUB>\","
738 <<
"\"N<SUB>refined</SUB>\","
739 <<
"\"N<SUB>unrefined</SUB>\","
740 <<
"\"norm of error\","
741 <<
"\"norm of solution\""
743 Trace_file <<
"ZONE T=\"" << time_stepper_pt()->type()
744 << time_stepper_pt()->nprev_values()
746 << time_pt()->dt() <<
"; ";
747 if (time_stepper_pt()->adaptive_flag())
749 Trace_file <<
"adaptive; target error " << Epsilon_t <<
" \"" << std::endl;
753 Trace_file <<
"(fixed timestep)\"" << std::endl;
763 template<
class ELEMENT>
774 cout <<
"=================================================" << std::endl;
775 cout <<
"Docing solution for t=" << time_pt()->time() << std::endl;
776 cout <<
"=================================================" << std::endl;
780 sprintf(filename,
"%s/soln%i.dat",Doc_info.directory().c_str(),
782 some_file.open(filename);
783 Bulk_mesh_pt->output(some_file,npts);
784 some_file <<
"TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
785 << time_pt()->time() <<
"\"";
786 some_file <<
"GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
787 some_file <<
"1" << std::endl;
788 some_file <<
"2" << std::endl;
789 some_file <<
" 0 0" << std::endl;
790 some_file << time_pt()->time()*20.0 <<
" 0" << std::endl;
793 some_file <<
"ZONE I=2,J=2" << std::endl;
794 some_file <<
"0.0 0.0 -1.2" << std::endl;
795 some_file <<
"1.3 0.0 -1.2" << std::endl;
796 some_file <<
"0.0 1.3 -1.2" << std::endl;
797 some_file <<
"1.3 1.3 -1.2" << std::endl;
798 some_file <<
"ZONE I=2,J=2" << std::endl;
799 some_file <<
"0.0 0.0 1.2" << std::endl;
800 some_file <<
"1.3 0.0 1.2" << std::endl;
801 some_file <<
"0.0 1.3 1.2" << std::endl;
802 some_file <<
"1.3 1.3 1.2" << std::endl;
809 sprintf(filename,
"%s/exact_soln%i.dat",Doc_info.directory().c_str(),
811 some_file.open(filename);
812 Bulk_mesh_pt->output_fct(some_file,npts,time_pt()->time(),
816 some_file <<
"ZONE I=2,J=2" << std::endl;
817 some_file <<
"0.0 0.0 -1.2" << std::endl;
818 some_file <<
"1.3 0.0 -1.2" << std::endl;
819 some_file <<
"0.0 1.3 -1.2" << std::endl;
820 some_file <<
"1.3 1.3 -1.2" << std::endl;
821 some_file <<
"ZONE I=2,J=2" << std::endl;
822 some_file <<
"0.0 0.0 1.2" << std::endl;
823 some_file <<
"1.3 0.0 1.2" << std::endl;
824 some_file <<
"0.0 1.3 1.2" << std::endl;
825 some_file <<
"1.3 1.3 1.2" << std::endl;
833 sprintf(filename,
"%s/error%i.dat",Doc_info.directory().c_str(),
835 some_file.open(filename);
836 Bulk_mesh_pt->compute_error(some_file,
846 cout <<
"error: " << error << std::endl;
847 cout <<
"norm : " << norm << std::endl << std::endl;
850 x[0]=Doc_node_pt->x(0);
851 x[1]=Doc_node_pt->x(1);
854 Vector<double > xi_wall(1);
855 Vector<double > r_wall(2);
857 Boundary_pt->position(xi_wall,r_wall);
858 Trace_file << time_pt()->time()
859 <<
" " << Doc_node_pt->value(0)
862 <<
" " << time_pt()->dt()
863 <<
" " << global_temporal_error_norm()
865 <<
" " << Bulk_mesh_pt->nelement()
866 <<
" " << Bulk_mesh_pt->nrefined()
867 <<
" " << Bulk_mesh_pt->nunrefined()
869 <<
" " << norm << std::endl;
873 sprintf(filename,
"%s/Wall%i.dat",Doc_info.directory().c_str(),
875 some_file.open(filename);
878 for (
unsigned iplot=0;iplot<nplot;iplot++)
880 xi_wall[0]=0.5*Pi*double(iplot)/double(nplot-1);
881 Boundary_pt->position(xi_wall,r_wall);
882 some_file << r_wall[0] <<
" " << r_wall[1] << std::endl;
887 sprintf(filename,
"%s/restart%i.dat",Doc_info.directory().c_str(),
889 some_file.open(filename);
905 template<
class ELEMENT>
908 Mesh*
const &surface_mesh_pt)
911 unsigned n_element = bulk_mesh_pt->nboundary_element(b);
914 for(
unsigned e=0;e<n_element;e++)
917 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
918 bulk_mesh_pt->boundary_element_pt(b,e));
921 int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
924 UnsteadyHeatFluxElement<ELEMENT>* flux_element_pt =
new
925 UnsteadyHeatFluxElement<ELEMENT>(bulk_elem_pt,face_index);
928 surface_mesh_pt->add_element_pt(flux_element_pt);
938 template<
class ELEMENT>
943 unsigned n_element = surface_mesh_pt->nelement();
946 for(
unsigned e=0;e<n_element;e++)
949 delete surface_mesh_pt->element_pt(e);
953 surface_mesh_pt->flush_element_and_node_storage();
962 template<
class ELEMENT>
965 double global_error = 0.0;
968 unsigned n_node = Bulk_mesh_pt->nnode();
971 for(
unsigned i=0;i<n_node;i++)
976 Bulk_mesh_pt->node_pt(i)->time_stepper_pt()->
977 temporal_error_in_value(Bulk_mesh_pt->node_pt(i),0);
980 global_error += error*error;
984 global_error /= double(n_node);
987 return sqrt(global_error);
995 template<
class ELEMENT>
1000 dump_file << Doc_info.number() <<
" # step number" << std::endl;
1003 dump_file << Next_dt <<
" # suggested next timestep" << std::endl;
1006 Problem::dump(dump_file);
1013 template<
class ELEMENT>
1017 string input_string;
1018 getline(restart_file,input_string,
'#');
1021 restart_file.ignore(80,
'\n');
1024 Doc_info.number()=unsigned(atof(input_string.c_str()));
1027 Doc_info.number()++;
1030 getline(restart_file,input_string,
'#');
1033 restart_file.ignore(80,
'\n');
1036 Next_dt=double(atof(input_string.c_str()));
1039 Problem::read(restart_file);
1058 CommandLineArgs::setup(argc,argv);
1078 unsigned max_adapt=10;
1087 cout <<
"Doing first timestep with dt: " << dt << std::endl;
1094 if (CommandLineArgs::Argc==2)
1102 problem.refine_uniformly();
1103 problem.refine_uniformly();
1111 while (problem.time_pt()->time()<t_max)
1115 problem.doubly_adaptive_unsteady_newton_solve(dt,problem.
epsilon_t(),
1117 cout <<
"Suggested new dt: " << dt_new << std::endl;
double B_hat
Amplitude of variation in y-half axis.
double A_hat
Amplitude of variation in x-half axis.
void position(const Vector< double > &xi, Vector< double > &r) const
Current position vector to material point at Lagrangian coordinate xi.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Parametrised position on object: r(xi). Evaluated at previous time level. t=0: current time; t>0: pre...
MyEllipse(const double &a, const double &b, const double &a_hat, const double &b_hat, const double &period, Time *time_pt)
Constructor: Pass half axes, amplitudes of their variation, period of oscillation and pointer to time...
double T
Period of oscillation.
Time * Time_pt
Pointer to time object.
virtual ~MyEllipse()
Destructor: Empty.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of prescribed flux elements.
void dump_it(ofstream &dump_file)
Dump problem data to allow for later restart.
Mesh * Surface_mesh_pt
Pointer to the "surface" mesh.
double Next_dt
Suggestion for next timestep (stored to allow it to be written to (or read from) restart file)
void set_initial_condition()
Set initial condition (incl previous timesteps) according to specified function.
double Epsilon_t
Target error for adaptive timestepping.
GeomObject * Boundary_pt
Pointer to GeomObject that specifies the domain bondary.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of prescribed flux elements.
RefineableQuarterCircleSectorMesh< ELEMENT > * bulk_mesh_pt()
Pointer to bulk mesh.
void create_flux_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &surface_mesh_pt)
Create UnsteadyHeat flux elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them t...
void doc_solution()
Doc the solution.
double & epsilon_t()
Target error for adaptive timestepping.
Node * Doc_node_pt
Pointer to central node (exists at all refinement levels) for doc.
RefineableUnsteadyHeatProblem(UnsteadyHeatEquations< 2 >::UnsteadyHeatSourceFctPt source_fct_pt)
Constructor: Pass pointer to source function and timestep.
double global_temporal_error_norm()
Global error norm for adaptive time-stepping.
ofstream Trace_file
Trace file.
void write_trace_file_header()
Write header for trace file.
~RefineableUnsteadyHeatProblem()
Destructor: Close trace file.
UnsteadyHeatEquations< 2 >::UnsteadyHeatSourceFctPt Source_fct_pt
Pointer to source function.
DocInfo Doc_info
Doc info object.
void actions_after_implicit_timestep()
Update the problem specs after timestep (empty)
void actions_before_newton_solve()
Update the problem specs before solve (empty)
void actions_before_implicit_timestep()
Update the problem specs before next timestep: Set Dirchlet boundary conditions from exact solution.
void delete_flux_elements(Mesh *const &surface_mesh_pt)
Delete UnsteadyHeat flux elements and wipe the surface mesh.
void actions_after_newton_solve()
Update the problem specs after solve (empty)
double & next_dt()
Suggestion for next timestep (stored to allow it to be written to (or read from) restart file)
void restart(ifstream &restart_file)
Read problem data for restart.
RefineableQuarterCircleSectorMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void get_exact_u(const double &time, const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
double Alpha
Parameter for steepness of step.
double Gamma
Parameter for timescale of step translation.
void get_exact_u(const double &time, const Vector< double > &x, double &u)
Exact solution as a scalar.
double Beta
Parameter for amplitude of step translation.
double step_position(const double &time)
Position of step (x-axis intercept)
void get_source(const double &time, const Vector< double > &x, double &source)
Source function to make it an exact solution.
void prescribed_flux_on_fixed_y_boundary(const double &time, const Vector< double > &x, double &flux)
Flux required by the exact solution on a boundary on which y is fixed.
double TanPhi
Parameter for angle of step.
int main(int argc, char *argv[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...