30 #include <feel/feel.hpp>
40 Feel::po::options_description
43 Feel::po::options_description stokesoptions(
"Stokes options" );
44 stokesoptions.add_options()
45 (
"penal", Feel::po::value<double>()->default_value( 0.5 ),
"penalisation parameter" )
46 (
"f", Feel::po::value<double>()->default_value( 0 ),
"forcing term" )
47 (
"mu", Feel::po::value<double>()->default_value( 1.0 ),
"reaction coefficient component" )
48 (
"hsize", Feel::po::value<double>()->default_value( 0.1 ),
"first h value to start convergence" )
49 (
"bctype", Feel::po::value<int>()->default_value( 0 ),
"0 = strong Dirichlet, 1 = weak Dirichlet" )
50 (
"bccoeff", Feel::po::value<double>()->default_value( 100.0 ),
"coeff for weak Dirichlet conditions" )
51 (
"geofile", Feel::po::value<std::string>()->default_value(
"backstep.geo" ),
"geometry file name" )
52 (
"export-matlab",
"export matrix and vectors in matlab" )
54 return stokesoptions.add( Feel::feel_options() ) ;
75 typedef double value_type;
77 typedef Backend<value_type> backend_type;
78 typedef boost::shared_ptr<backend_type> backend_ptrtype;
81 typedef Simplex<FEELPP_NS_DIM> convex_type;
82 typedef Mesh<convex_type> mesh_type;
83 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
87 typedef Lagrange<FEELPP_NS_ORDER+1, Vectorial> basis_u_type;
88 typedef Lagrange<FEELPP_NS_ORDER, Scalar> basis_p_type;
89 typedef Lagrange<0, Scalar> basis_l_type;
91 #if defined( FEELPP_USE_LM )
92 typedef bases<basis_u_type,basis_p_type, basis_l_type> basis_type;
94 typedef bases<basis_u_type,basis_p_type> basis_type;
100 typedef FunctionSpace<mesh_type, basis_type> space_type;
101 typedef boost::shared_ptr<space_type> space_ptrtype;
106 typedef space_type::element_type element_type;
109 typedef BDF<space_type> bdf_type;
110 typedef boost::shared_ptr<bdf_type> bdf_ptrtype;
113 typedef Exporter<mesh_type> export_type;
135 void exportResults(
double time, element_type& u );
139 backend_ptrtype M_backend;
148 boost::shared_ptr<export_type> exporter;
152 NavierStokes::NavierStokes()
155 M_backend( backend_type::build( this->vm() ) ),
156 meshSize( this->vm()[
"hsize"].as<double>() ),
157 mu( this->vm()[
"mu"].as<value_type>() ),
158 penalbc( this->vm()[
"bccoeff"].as<value_type>() ),
159 exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) )
167 Environment::changeRepository( boost::format(
"doc/manual/ns/%1%/%2%/P%3%P%4%/h_%5%/" )
168 % this->about().appName()
169 % convex_type::name()
170 % basis_u_type::nOrder % basis_p_type::nOrder
171 % this->vm()[
"hsize"].as<double>() );
174 mesh = createGMSHMesh( _mesh=
new mesh_type,
175 _desc=geo( _filename=this->vm()[
"geofile"].as<std::string>()
180 Xh = space_type::New( mesh );
188 auto U = Xh->element(
"(u,p)" );
189 auto V = Xh->element(
"(u,q)" );
190 auto u = U.element<0>(
"u" );
191 auto v = V.element<0>(
"u" );
192 auto p = U.element<1>(
"p" );
193 auto q = V.element<1>(
"p" );
194 #if defined( FEELPP_USE_LM )
195 auto lambda = U.element<2>();
196 auto nu = V.element<2>();
200 LOG(INFO) <<
"[dof] number of dof: " << Xh->nDof() <<
"\n";
201 LOG(INFO) <<
"[dof] number of dof/proc: " << Xh->nLocalDof() <<
"\n";
202 LOG(INFO) <<
"[dof] number of dof(U): " << Xh->functionSpace<0>()->nDof() <<
"\n";
203 LOG(INFO) <<
"[dof] number of dof/proc(U): " << Xh->functionSpace<0>()->nLocalDof() <<
"\n";
204 LOG(INFO) <<
"[dof] number of dof(P): " << Xh->functionSpace<1>()->nDof() <<
"\n";
205 LOG(INFO) <<
"[dof] number of dof/proc(P): " << Xh->functionSpace<1>()->nLocalDof() <<
"\n";
207 LOG(INFO) <<
"Data Summary:\n";
208 LOG(INFO) <<
" hsize = " << meshSize <<
"\n";
209 LOG(INFO) <<
" export = " << this->vm().count(
"export" ) <<
"\n";
210 LOG(INFO) <<
" mu = " << mu <<
"\n";
211 LOG(INFO) <<
" bccoeff = " << penalbc <<
"\n";
217 auto deft = gradt( u )+trans(gradt(u));
218 auto def = grad( v )+trans(grad(v));
223 auto SigmaNt = -idt( p )*N()+mu*deft*N();
226 auto SigmaN = -id( p )*N()+mu*def*N();
229 auto F = M_backend->newVector( Xh );
230 auto D = M_backend->newMatrix( Xh, Xh );
233 auto ns_rhs = form1( _test=Xh, _vector=F );
236 LOG(INFO) <<
"[navier-stokes] vector local assembly done\n";
239 auto bdfns=bdf(_space=Xh);
245 auto navierstokes = form2( _test=Xh, _trial=Xh, _matrix=D );
247 navierstokes +=
integrate(
elements( mesh ), mu*inner( deft,def )+ trans(idt( u ))*
id( v )*bdfns->polyDerivCoefficient( 0 ) );
248 LOG(INFO) <<
"mu*inner(deft,def)+(bdf(u),v): " << chrono.elapsed() <<
"\n";
250 navierstokes +=
integrate(
elements( mesh ), - div( v )*idt( p ) + divt( u )*
id( q ) );
251 LOG(INFO) <<
"(u,p): " << chrono.elapsed() <<
"\n";
253 #if defined( FEELPP_USE_LM )
254 navierstokes +=
integrate(
elements( mesh ),
id( q )*idt( lambda ) + idt( p )*
id( nu ) );
255 LOG(INFO) <<
"(lambda,p): " << chrono.elapsed() <<
"\n";
259 std::for_each( inflow_conditions.begin(), inflow_conditions.end(),
260 [&]( BoundaryCondition
const& bc )
263 ns_rhs +=
integrate(
markedfaces( mesh, bc.marker() ), inner( idf(&bc,BoundaryCondition::operator()),-SigmaN+penalbc*id( v )/hFace() ) );
269 std::for_each( wall_conditions.begin(), wall_conditions.end(),
270 [&]( BoundaryCondition
const& bc )
276 std::for_each( outflow_conditions.begin(), outflow_conditions.end(),
277 [&]( BoundaryCondition
const& bc )
279 ns_rhs +=
integrate(
markedfaces( mesh, bc.marker() ), inner( idf(&bc,BoundaryCondition::operator()),N() ) );
282 LOG(INFO) <<
"bc: " << chrono.elapsed() <<
"\n";
285 u =
vf::project( _space=Xh->functionSpace<0>(), _expr=cst(0.) );
286 p =
vf::project( _space=Xh->functionSpace<1>(), _expr=cst(0.) );
288 M_bdf->initialize( U );
290 for( bdfns->start(); bdfns->isFinished(); bdfns->next() )
293 auto bdf_poly = bdfns->polyDeriv();
294 form1( _test=Xh, _vector=Ft ) =
297 form1( _test=Xh, _vector=Ft ) +=
303 backend()->solve( _matrix=D, _solution=U, _rhs=Ft );
305 this->exportResults( bdfns->time(), U );
317 Navierstokes::exportResults(
double time, element_type& U )
319 auto u = U.element<0>();
320 auto p = U.element<1>();
322 auto v = V.element<0>();
323 auto q = V.element<1>();
324 #if defined( FEELPP_USE_LM )
325 auto lambda = U.element<2>();
326 auto nu = V.element<2>();
327 LOG(INFO) <<
"value of the Lagrange multiplier lambda= " << lambda( 0 ) <<
"\n";
331 double div_u_error_L2 = normL2( _range=
elements( u.mesh() ), _expr=divv( u ) );
332 LOG(INFO) <<
"[navierstokes] ||div(u)||_2=" << div_u_error_L2 <<
"\n";
334 if ( exporter->doExport() )
336 exporter->step( time )->setMesh( U.functionSpace()->mesh() );
337 exporter->step( time )->addRegions();
338 exporter->step( time )->add( {
"u",
"p",
"l"}, U );
virtual void run(const double *X, unsigned long P, double *Y, unsigned long N)
Definition: simget.hpp:202
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::marker_face_const_iterator, typename MeshTraits< MeshType >::marker_face_const_iterator > markedfaces(MeshType const &mesh)
Definition: filters.hpp:969
Poly::value_type integrate(Polynomial< Poly, PsetType > const &p)
integrate a polynomial over a convex
Definition: operations.hpp:51
elements_type const & elements() const
Definition: elements.hpp:355
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
double meshSize() const
return the mesh size
Definition: simget.hpp:137
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::location_face_const_iterator, typename MeshTraits< MeshType >::location_face_const_iterator > boundaryfaces(MeshType const &mesh)
Definition: filters.hpp:1158
Simget()
Definition: simget.cpp:36