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() ) );
265 navierstokes +=integrate( boundaryfaces( mesh ), -inner( SigmaNt,
id( v ) ) );
266 navierstokes +=integrate( boundaryfaces( mesh ), -inner( SigmaN,idt( u ) ) );
267 navierstokes +=integrate( boundaryfaces( mesh ), +penalbc*inner( idt( u ),
id( v ) )/hFace() );
269 std::for_each( wall_conditions.begin(), wall_conditions.end(),
270 [&]( BoundaryCondition
const& bc )
272 navierstokes +=integrate( boundaryfaces( mesh ), -inner( SigmaNt,
id( v ) ) );
273 navierstokes +=integrate( boundaryfaces( mesh ), -inner( SigmaN,idt( u ) ) );
274 navierstokes +=integrate( boundaryfaces( mesh ), +penalbc*inner( idt( u ),
id( v ) )/hFace() );
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 ) =
295 integrate( _range=elements(mesh), _expr=trans(idv( bdf_poly ))*
id( v ) );
297 form1( _test=Xh, _vector=Ft ) +=
298 integrate( _range=elements(mesh), _expr=trans(gradv(u)*idv( u ))*
id(v) );
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 );
Feel::po::options_description makeOptions()
Definition: ns.hpp:41