29 #ifndef FEELPP_DRIVENCAVITY_IMPL_HPP_H
30 #define FEELPP_DRIVENCAVITY_IMPL_HPP_H 1
37 DrivenCavity<Dim>::DrivenCavity( )
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() ) )
48 void DrivenCavity<Dim>::init()
50 if ( this->vm().count(
"help" ) )
52 if ( Environment::worldComm().isMasterRank() )
53 std::cout << this->optionsDescription() <<
"\n";
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";
60 Vh = space_type::New( mesh );
65 DrivenCavity<Dim>::Jacobian(
const vector_ptrtype& X, sparse_matrix_ptrtype& J)
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" );
74 auto lambda = U.template element<2>();
75 auto nu = V.template element<2>();
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) );
84 a += integrate( elements( mesh ), trans(
id(v))*gradv(u)*idt(u));
85 a += integrate( elements( mesh ), trans(
id(v))*gradt(u)*idv(u));
88 a += integrate(elements(mesh),
id(q)*idt(lambda)+idt(p)*
id(nu));
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() );
102 DrivenCavity<Dim>::Residual(
const vector_ptrtype& X, vector_ptrtype& R)
104 auto U = Vh->element(
"(u,p)" );
105 auto V = Vh->element(
"(v,q)" );
106 auto u = U.template element<0>(
"u" );
108 auto v = V.template element<0>(
"u" );
109 auto p = U.template element<1>(
"p" );
110 auto q = V.template element<1>(
"p" );
112 auto lambda = U.template element<2>();
113 auto nu = V.template element<2>();
117 auto u_exact=vf::project(Vh->template functionSpace<0>(), markedfaces(mesh,
"wall2"), uex );
121 auto r = form1( _test=Vh, _vector=R );
123 r = integrate( elements( mesh ), trans(gradv( u )*idv(u))*
id(v));
124 r += integrate( elements( mesh ), inner(gradv( u ),grad( v ))/Re );
125 r += integrate( elements( mesh ),-idv(p)*div(v) +
id(q)*divv(u));
127 r += integrate ( elements( mesh ), +
id( q )*idv( lambda )+idv( p )*
id( nu ) );
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() );
140 void DrivenCavity<Dim>::exportResults( element_type
const& U )
143 auto u_exact=vf::project(Vh->template functionSpace<0>(), markedfaces(mesh,
"wall2"), uex );
145 if ( exporter->doExport() )
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);
157 void DrivenCavity<Dim>::run()
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" );
168 auto lambda = U.template element<2>();
169 auto nu = V.template element<2>();
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));
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 );
186 if ( option(_name=
"continuation" ).
template as<bool>() )
188 double ReTarget = Re;
189 int N = std::ceil( std::log( Re ) );
190 for(
int i = 0; i <= N; ++i )
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";
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";
209 this->exportResults( U );