Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Files Functions Variables Typedefs Pages
drivencavity_impl.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: 2013-11-22
7 
8  Copyright (C) 2013 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 */
29 #ifndef FEELPP_DRIVENCAVITY_IMPL_HPP_H
30 #define FEELPP_DRIVENCAVITY_IMPL_HPP_H 1
31 
32 #include "drivencavity.hpp"
33 
34 namespace Feel
35 {
36 template<int Dim>
37 DrivenCavity<Dim>::DrivenCavity( )
38  :
39  super( ),
40  Re( this->vm()["Re"].template as<value_type>() ),
41  penalbc( this->vm()["bccoeff"].template as<value_type>() ),
42  exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) )
43 {
44 
45 }
46 
47 template<int Dim>
48 void DrivenCavity<Dim>::init()
49 {
50  if ( this->vm().count( "help" ) )
51  {
52  if ( Environment::worldComm().isMasterRank() )
53  std::cout << this->optionsDescription() << "\n";
54  return;
55  }
56  mesh = loadMesh(_mesh=new Mesh<Simplex<Dim>>);
57  if ( Environment::worldComm().isMasterRank() )
58  std::cout << "number of elements of " << Dim << "D: " << mesh->numElements() << "\n";
59 
60  Vh = space_type::New( mesh );
61 }
62 
63 template<int Dim>
64 void
65 DrivenCavity<Dim>::Jacobian(const vector_ptrtype& X, sparse_matrix_ptrtype& J)
66 {
67  auto U = Vh->element( "(u,p)" );
68  auto V = Vh->element( "(v,q)" );
69  auto u = U.template element<0>( "u" );
70  auto v = V.template element<0>( "u" );
71  auto p = U.template element<1>( "p" );
72  auto q = V.template element<1>( "p" );
73  //#if defined( FEELPP_USE_LM )
74  auto lambda = U.template element<2>();
75  auto nu = V.template element<2>();
76 
77  //#endif
78 
79  if (!J) J = backend(_name="newtonns")->newMatrix( Vh, Vh );
80  auto a = form2( _test=Vh, _trial=Vh, _matrix=J );
81  a = integrate( elements( mesh ), inner(gradt( u ),grad( v ))/Re );
82  a += integrate( elements( mesh ), id(q)*divt(u) -idt(p)*div(v) );
83  // Convective terms
84  a += integrate( elements( mesh ), trans(id(v))*gradv(u)*idt(u));
85  a += integrate( elements( mesh ), trans(id(v))*gradt(u)*idv(u));
86 
87  //#if defined( FEELPP_USE_LM )
88  a += integrate(elements(mesh), id(q)*idt(lambda)+idt(p)*id(nu));
89  //#elif
90  //a += integrate(elements(mesh), idt(p)*id(nu));
91 
92  //Weak Dirichlet conditions
93  a += integrate( boundaryfaces( mesh ),-trans( -idt(p)*N()+gradt(u)*N()/Re )*id( v ));//
94  a += integrate( boundaryfaces( mesh ),-trans( -id(p)*N()+grad(u)*N()/Re )*idt( u ));//
95  a += integrate( boundaryfaces( mesh ), +penalbc*inner( idt( u ),id( v ) )/hFace() );
96 
97 
98 }
99 
100 template<int Dim>
101 void
102 DrivenCavity<Dim>::Residual(const vector_ptrtype& X, vector_ptrtype& R)
103 {
104  auto U = Vh->element( "(u,p)" );
105  auto V = Vh->element( "(v,q)" );
106  auto u = U.template element<0>( "u" );
107  //auto u_exact = U.template element<0>( "u_exact" );
108  auto v = V.template element<0>( "u" );
109  auto p = U.template element<1>( "p" );
110  auto q = V.template element<1>( "p" );
111  //#if defined( FEELPP_USE_LM )
112  auto lambda = U.template element<2>();
113  auto nu = V.template element<2>();
114  //#endif
115 
116  auto uex=unitX();
117  auto u_exact=vf::project(Vh->template functionSpace<0>(), markedfaces(mesh, "wall2"), uex );
118 
119 
120  U=*X;
121  auto r = form1( _test=Vh, _vector=R );
122  //r += integrate( elements( mesh ),-inner( f,id( v ) ) );
123  r = integrate( elements( mesh ), trans(gradv( u )*idv(u))*id(v));//convective term
124  r += integrate( elements( mesh ), inner(gradv( u ),grad( v ))/Re );
125  r += integrate( elements( mesh ),-idv(p)*div(v) + id(q)*divv(u));
126  //#if defined( FEELPP_USE_LM )
127  r += integrate ( elements( mesh ), +id( q )*idv( lambda )+idv( p )*id( nu ) );
128  //#endif
129 
130 
131  //Weak Dirichlet
132  auto SigmaNv = ( -idv( p )*N() + gradv( u )*N()/Re );
133  auto SigmaN = ( -id( q )*N() + grad( v )*N()/Re );
134  r +=integrate ( boundaryfaces(mesh), - trans( SigmaNv )*id( v ) - trans( SigmaN )*( idv( u ) -idv(u_exact) ) + penalbc*trans( idv( u ) - idv(u_exact) )*id( v )/hFace() );
135 
136 
137 }
138 
139 template<int Dim>
140 void DrivenCavity<Dim>::exportResults( element_type const& U )
141 {
142  auto uex=unitX();
143  auto u_exact=vf::project(Vh->template functionSpace<0>(), markedfaces(mesh, "wall2"), uex );
144 
145  if ( exporter->doExport() )
146  {
147  exporter->step( 0 )->setMesh( U.functionSpace()->mesh() );
148  exporter->step( 0 )->addRegions();
149  exporter->step( 0 )->add( "u", U.template element<0>() );
150  exporter->step( 0 )->add( "p", U.template element<1>() );
151  exporter->step( 0 )->add( "uex", u_exact);
152  exporter->save();
153  }
154 
155 }
156 template<int Dim>
157 void DrivenCavity<Dim>::run()
158 {
159  this->init();
160 
161  auto U = Vh->element( "(u,p)" );
162  auto V = Vh->element( "(v,q)" );
163  auto u = U.template element<0>( "u" );
164  auto v = V.template element<0>( "u" );
165  auto p = U.template element<1>( "p" );
166  auto q = V.template element<1>( "p" );
167  //#if defined( FEELPP_USE_LM )
168  auto lambda = U.template element<2>();
169  auto nu = V.template element<2>();
170  //#endif
171 
172 
173 
174  u=vf::project(Vh->template functionSpace<0>(), elements(mesh), zero<Dim,1>());
175  p=vf::project(Vh->template functionSpace<1>(), elements(mesh), constant(0.0));
176 
177  if ( Environment::worldComm().isMasterRank() )
178  std::cout << "Initializing residual " << "\n";
179  backend(_name="newtonns")->nlSolver()->residual = boost::bind( &DrivenCavity::Residual,
180  boost::ref( *this ), _1, _2 );
181  if ( Environment::worldComm().isMasterRank() )
182  std::cout << "Initializing the Jacobian matrix" << "\n";
183  backend(_name="newtonns")->nlSolver()->jacobian = boost::bind( &DrivenCavity::Jacobian,
184  boost::ref( *this ), _1, _2 );
185 
186  if ( option(_name="continuation" ).template as<bool>() )
187  {
188  double ReTarget = Re;
189  int N = std::ceil( std::log( Re ) );
190  for( int i = 0; i <= N; ++i )
191  {
192  Re = std::exp( std::log(1)+i*(std::log(ReTarget)-std::log(1))/double(N));
193  if ( Environment::worldComm().isMasterRank() )
194  std::cout << "Start solving for Reynolds = " << Re << "\n";
195  backend(_name="newtonns")->nlSolve( _solution=U );
196  if ( Environment::worldComm().isMasterRank() )
197  std::cout << "Done!" << "\n";
198  }
199  }
200  else
201  {
202  if ( Environment::worldComm().isMasterRank() )
203  std::cout << "Start solving for Reynolds = " << Re << "\n";
204  backend(_name="newtonns")->nlSolve( _solution=U );
205  if ( Environment::worldComm().isMasterRank() )
206  std::cout << "Done!" << "\n";
207  }
208 
209  this->exportResults( U );
210 };
211 } // namespace Feel
212 #endif /* FEELPP_DRIVENCAVITY_IMPL_HPP_H */