32 #include "advection_diffusion.h"
33 #include "navier_stokes.h"
34 #include "multi_physics.h"
38 #include "meshes/rectangular_quadmesh.h"
41 using namespace oomph;
79 template<
class NST_ELEMENT,
class AD_ELEMENT>
92 void actions_before_newton_solve();
102 return dynamic_cast<RefineableRectangularQuadMesh<NST_ELEMENT>*
>
111 return dynamic_cast<RefineableRectangularQuadMesh<AD_ELEMENT>*
>
123 RefineableNavierStokesEquations<2>::
124 unpin_all_pressure_dofs(nst_mesh_pt()->element_pt());
128 fix_pressure(0,0,0.0);
131 Multi_domain_functions::
132 setup_multi_domain_interactions<NST_ELEMENT,AD_ELEMENT>
133 (
this,nst_mesh_pt(),adv_diff_mesh_pt());
140 const double &pvalue)
143 if (nst_mesh_pt()->nelement()>0)
145 dynamic_cast<NST_ELEMENT*
>(nst_mesh_pt()->element_pt(e))->
146 fix_pressure(pdof,pvalue);
185 template<
class NST_ELEMENT,
class AD_ELEMENT>
190 Doc_info.set_directory(
"RESLT");
206 new RefineableRectangularQuadMesh<NST_ELEMENT>(n_x,n_y,l_x,l_y);
208 new RefineableRectangularQuadMesh<AD_ELEMENT>(n_x,n_y,l_x,l_y);
211 Nst_mesh_pt->spatial_error_estimator_pt()=
new Z2ErrorEstimator;
212 Adv_diff_mesh_pt->spatial_error_estimator_pt()=
new Z2ErrorEstimator;
215 Nst_mesh_pt->max_permitted_error()=0.5e-3;
216 Nst_mesh_pt->min_permitted_error()=0.5e-4;
217 Adv_diff_mesh_pt->max_permitted_error()=0.5e-3;
218 Adv_diff_mesh_pt->min_permitted_error()=0.5e-4;
226 unsigned num_bound = nst_mesh_pt()->nboundary();
227 for(
unsigned ibound=0;ibound<num_bound;ibound++)
233 unsigned num_nod= nst_mesh_pt()->nboundary_node(ibound);
234 for (
unsigned inod=0;inod<num_nod;inod++)
240 if((ibound==1) || (ibound==3))
246 val_max=nst_mesh_pt()->boundary_node_pt(ibound,inod)->nvalue();
250 for(
unsigned j=0;j<val_max;j++)
252 nst_mesh_pt()->boundary_node_pt(ibound,inod)->pin(j);
259 fix_pressure(0,0,0.0);
262 num_bound = adv_diff_mesh_pt()->nboundary();
263 for(
unsigned ibound=0;ibound<num_bound;ibound++)
269 unsigned num_nod= adv_diff_mesh_pt()->nboundary_node(ibound);
270 for (
unsigned inod=0;inod<num_nod;inod++)
275 if ((ibound==1) || (ibound==3))
281 val_max=adv_diff_mesh_pt()->boundary_node_pt(ibound,inod)->nvalue();
283 for(
unsigned j=0;j<val_max;j++)
285 adv_diff_mesh_pt()->boundary_node_pt(ibound,inod)->pin(j);
296 unsigned n_nst_element = nst_mesh_pt()->nelement();
297 for(
unsigned i=0;i<n_nst_element;i++)
300 NST_ELEMENT *el_pt =
dynamic_cast<NST_ELEMENT*
>
301 (nst_mesh_pt()->element_pt(i));
319 el_pt->ignore_external_geometric_data();
322 unsigned n_ad_element = adv_diff_mesh_pt()->nelement();
323 for(
unsigned i=0;i<n_ad_element;i++)
326 AD_ELEMENT *el_pt =
dynamic_cast<AD_ELEMENT*
>
327 (adv_diff_mesh_pt()->element_pt(i));
339 el_pt->ignore_external_geometric_data();
344 add_sub_mesh(Nst_mesh_pt);
345 add_sub_mesh(Adv_diff_mesh_pt);
349 Multi_domain_functions::
350 setup_multi_domain_interactions<NST_ELEMENT,AD_ELEMENT>
351 (
this,nst_mesh_pt(),adv_diff_mesh_pt());
354 cout <<
"Number of equations: " << assign_eqn_numbers() << endl;
363 template<
class NST_ELEMENT,
class AD_ELEMENT>
367 unsigned num_bound = nst_mesh_pt()->nboundary();
368 for(
unsigned ibound=0;ibound<num_bound;ibound++)
371 unsigned num_nod=nst_mesh_pt()->nboundary_node(ibound);
372 for(
unsigned inod=0;inod<num_nod;inod++)
375 Node* nod_pt=nst_mesh_pt()->boundary_node_pt(ibound,inod);
380 if((ibound==1) || (ibound==3)) {vel_max = 1;}
382 for(
unsigned j=0;j<vel_max;j++) {nod_pt->set_value(j,0.0);}
391 double x = nod_pt->x(0);
394 double value = sin(2.0*3.141592654*x/3.0);
395 nod_pt->set_value(1,value);
403 num_bound=adv_diff_mesh_pt()->nboundary();
404 for(
unsigned ibound=0;ibound<num_bound;ibound++)
407 unsigned num_nod=adv_diff_mesh_pt()->nboundary_node(ibound);
408 for(
unsigned inod=0;inod<num_nod;inod++)
411 Node* nod_pt=adv_diff_mesh_pt()->boundary_node_pt(ibound,inod);
415 if(ibound==2) {nod_pt->set_value(0,-0.5);}
419 if(ibound==0) {nod_pt->set_value(0,0.5);}
431 template<
class NST_ELEMENT,
class AD_ELEMENT>
442 sprintf(filename,
"%s/fluid_soln%i.dat",Doc_info.directory().c_str(),
444 some_file.open(filename);
445 nst_mesh_pt()->output(some_file,npts);
449 sprintf(filename,
"%s/temperature_soln%i.dat",Doc_info.directory().c_str(),
451 some_file.open(filename);
452 adv_diff_mesh_pt()->output(some_file,npts);
475 RefineableQCrouzeixRaviartElement<2>,
476 RefineableQAdvectionDiffusionElement<2,3> > ,
478 RefineableQAdvectionDiffusionElement<2,3>,
479 RefineableQCrouzeixRaviartElement<2> >
484 problem.enable_imperfection();
487 problem.newton_solve(2);
490 problem.doc_solution();
497 problem.disable_imperfection();
498 problem.newton_solve(2);
499 problem.doc_solution();
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
void actions_after_newton_solve()
Update the problem after solve (empty)
void fix_pressure(const unsigned &e, const unsigned &pdof, const double &pvalue)
Fix pressure in element e at pressure dof pdof and set to pvalue.
~RefineableConvectionProblem()
Destructor. Empty.
void doc_solution()
Doc the solution.
RefineableRectangularQuadMesh< AD_ELEMENT > * adv_diff_mesh_pt()
Access function to the AD mesh. Casts the pointer to the base Mesh object to the actual mesh type.
void actions_before_newton_solve()
Update the problem specs before solve:
void disable_imperfection()
Set th boundary condition on the upper wall to be unperturbed.
void actions_before_adapt()
Actions before adapt:(empty)
RefineableConvectionProblem()
Constructor.
void actions_after_adapt()
Actions after adaptation, reset all sources, then re-pin a single pressure degree of freedom.
RefineableRectangularQuadMesh< NST_ELEMENT > * nst_mesh_pt()
Access function to the NST mesh. Casts the pointer to the base Mesh object to the actual mesh type.
void enable_imperfection()
Set the boundary condition on the upper wall to be perturbed slightly to force the solution into the ...
RefineableRectangularQuadMesh< NST_ELEMENT > * Nst_mesh_pt
Navier Stokes mesh.
RefineableRectangularQuadMesh< AD_ELEMENT > * Adv_diff_mesh_pt
Advection diffusion mesh.
DocInfo Doc_info
DocInfo object.
bool Imperfect
Is the boundary condition imperfect or not.
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
int main()
Driver code for 2D Boussinesq convection problem with adaptivity.
Namespace for the physical parameters in the problem.
Vector< double > Direction_of_gravity(2)
Gravity vector.
double Rayleigh
Rayleigh number, set to be greater than the threshold for linear instability.
double Inverse_Prandtl
1/Prandtl number
double Peclet
Peclet number (identically one from our non-dimensionalisation)