Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ns.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*-
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <prudhomme@unistra.fr>
6  Date: 2012-09-30
7 
8  Copyright (C) 2012 Université de Strasbourg
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 2.1 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
30 #include <feel/feel.hpp>
31 
39 inline
40 Feel::po::options_description
41 makeOptions()
42 {
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" )
53  ;
54  return stokesoptions.add( Feel::feel_options() ) ;
55 }
56 
57 
58 
59 
60 namespace Feel
61 {
67 class NavierStokes
68  :
69 public Simget
70 {
71  typedef Simget super;
72 public:
73 
74 
75  typedef double value_type;
76 
77  typedef Backend<value_type> backend_type;
78  typedef boost::shared_ptr<backend_type> backend_ptrtype;
79 
80  /*mesh*/
81  typedef Simplex<FEELPP_NS_DIM> convex_type;
82  typedef Mesh<convex_type> mesh_type;
83  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
84 
85  /*basis*/
86  //# marker1 #,
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;
90  // use lagrange multipliers to ensure zero mean pressure
91 #if defined( FEELPP_USE_LM )
92  typedef bases<basis_u_type,basis_p_type, basis_l_type> basis_type;
93 #else
94  typedef bases<basis_u_type,basis_p_type> basis_type;
95 #endif
96  //# endmarker1 #
97 
98  /*space*/
99  //# marker2 #
100  typedef FunctionSpace<mesh_type, basis_type> space_type;
101  typedef boost::shared_ptr<space_type> space_ptrtype;
102  //# endmarker2 #
103 
104  /* functions */
105  //# marker3 #
106  typedef space_type::element_type element_type;
107  //# endmarker3 #
108 
109  typedef BDF<space_type> bdf_type;
110  typedef boost::shared_ptr<bdf_type> bdf_ptrtype;
111 
112  /* export */
113  typedef Exporter<mesh_type> export_type;
114 
115  FEELPP_DONT_INLINE
116  NavierStokes();
117 
118  // init mesh and space
119  FEELPP_DONT_INLINE
120  void init();
121 
125  FEELPP_DONT_INLINE
126  void run();
127 
128 private:
129 
130 
134  FEELPP_DONT_INLINE
135  void exportResults( double time, element_type& u );
136 
137 private:
138 
139  backend_ptrtype M_backend;
140  double meshSize;
141 
142  double mu;
143  double penalbc;
144 
145  mesh_ptrtype mesh;
146  space_ptrtype Xh;
147 
148  boost::shared_ptr<export_type> exporter;
149 }; // NavierStokes
150 
151 
152 NavierStokes::NavierStokes()
153  :
154  super(),
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() ) )
160 {
161 
162 }
163 
164 void
165 NavierStokes::init()
166 {
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>() );
172 
173 
174  mesh = createGMSHMesh( _mesh=new mesh_type,
175  _desc=geo( _filename=this->vm()["geofile"].as<std::string>()
176  _dim=FEELPP_NS_DIM,
177  _h=meshSize ) );
178 
179 
180  Xh = space_type::New( mesh );
181 
182 }
183 void
184 NavierStokes::run()
185 {
186  this->init();
187 
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>();
197 #endif
198  //# endmarker4 #
199 
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";
206 
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";
212 
213 
214 
215 
216  //# marker5 #
217  auto deft = gradt( u )+trans(gradt(u));
218  auto def = grad( v )+trans(grad(v));
219  //# endmarker5 #
220 
221  //# marker6 #
222  // total stress tensor (trial)
223  auto SigmaNt = -idt( p )*N()+mu*deft*N();
224 
225  // total stress tensor (test)
226  auto SigmaN = -id( p )*N()+mu*def*N();
227  //# endmarker6 #
228 
229  auto F = M_backend->newVector( Xh );
230  auto D = M_backend->newMatrix( Xh, Xh );
231 
232  // right hand side
233  auto ns_rhs = form1( _test=Xh, _vector=F );
234 
235 
236  LOG(INFO) << "[navier-stokes] vector local assembly done\n";
237 
238  // construction of the BDF
239  auto bdfns=bdf(_space=Xh);
240 
241  /*
242  * Construction of the left hand side
243  */
244 
245  auto navierstokes = form2( _test=Xh, _trial=Xh, _matrix=D );
246  mpi::timer chrono;
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";
249  chrono.restart();
250  navierstokes +=integrate( elements( mesh ), - div( v )*idt( p ) + divt( u )*id( q ) );
251  LOG(INFO) << "(u,p): " << chrono.elapsed() << "\n";
252  chrono.restart();
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";
256  chrono.restart();
257 #endif
258 
259  std::for_each( inflow_conditions.begin(), inflow_conditions.end(),
260  [&]( BoundaryCondition const& bc )
261  {
262  // right hand side
263  ns_rhs += integrate( markedfaces( mesh, bc.marker() ), inner( idf(&bc,BoundaryCondition::operator()),-SigmaN+penalbc*id( v )/hFace() ) );
264 
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() );
268  });
269  std::for_each( wall_conditions.begin(), wall_conditions.end(),
270  [&]( BoundaryCondition const& bc )
271  {
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() );
275  });
276  std::for_each( outflow_conditions.begin(), outflow_conditions.end(),
277  [&]( BoundaryCondition const& bc )
278  {
279  ns_rhs += integrate( markedfaces( mesh, bc.marker() ), inner( idf(&bc,BoundaryCondition::operator()),N() ) );
280  });
281 
282  LOG(INFO) << "bc: " << chrono.elapsed() << "\n";
283  chrono.restart();
284 
285  u = vf::project( _space=Xh->functionSpace<0>(), _expr=cst(0.) );
286  p = vf::project( _space=Xh->functionSpace<1>(), _expr=cst(0.) );
287 
288  M_bdf->initialize( U );
289 
290  for( bdfns->start(); bdfns->isFinished(); bdfns->next() )
291  {
292  // add time dependent terms
293  auto bdf_poly = bdfns->polyDeriv();
294  form1( _test=Xh, _vector=Ft ) =
295  integrate( _range=elements(mesh), _expr=trans(idv( bdf_poly ))*id( v ) );
296  // add convective terms
297  form1( _test=Xh, _vector=Ft ) +=
298  integrate( _range=elements(mesh), _expr=trans(gradv(u)*idv( u ))*id(v) );
299  // add contrib from time independent terms
300  Ft->add( 1., F );
301 
302  // add time stepping terms from BDF to right hand side
303  backend()->solve( _matrix=D, _solution=U, _rhs=Ft );
304 
305  this->exportResults( bdfns->time(), U );
306  }
307 
308 
309 
310 
311 
312 } // NavierNavierstokes::run
313 
314 
315 
316 void
317 Navierstokes::exportResults( double time, element_type& U )
318 {
319  auto u = U.element<0>();
320  auto p = U.element<1>();
321 
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";
328 
329 #endif
330 
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";
333 
334  if ( exporter->doExport() )
335  {
336  exporter->step( time )->setMesh( U.functionSpace()->mesh() );
337  exporter->step( time )->addRegions();
338  exporter->step( time )->add( {"u","p","l"}, U );
339  exporter->save();
340  }
341 
342 } // NavierStokes::export
343 } // Feel
344 
345 
346 
347 

Generated on Sun Oct 20 2013 08:25:02 for Feel++ by doxygen 1.8.4