13 #include <feel/feeldiscr/functionspace.hpp>
17 #include <feel/feelfilters/exporter.hpp>
29 #include <Eigen/Dense>
30 #include <Eigen/Eigenvalues>
35 #include<boost/range/algorithm/max_element.hpp>
36 #include<boost/filesystem.hpp>
40 using namespace Feel::vf;
44 template<
int PolynomialOrder>
class NIRBTEST
54 typedef double value_type;
59 typedef boost::shared_ptr<backend_type> backend_ptrtype;
63 typedef typename backend_type::vector_ptrtype vector_ptrtype;
67 typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
75 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
79 typedef bases<Lagrange<PolynomialOrder,Scalar> > basis_type;
84 typedef boost::shared_ptr<space_type> space_ptrtype;
85 typedef typename space_type::element_type element_type;
90 typedef boost::shared_ptr<export_type> export_ptrtype;
91 typedef std::vector<element_type> vector_of_element_type;
100 CoarseMeshSize( this->vm()[
"hcoarsesize"].template as<double>() ),
101 FineMeshSize( this->vm()[
"hfinsize"].template as<double>() ),
102 ReadingMeshes( this->vm()[
"ReadingMeshes"].template as<int>() ),
103 NbSnapshot( this->vm()[
"NbSnapshot"].template as<int>() ),
104 sizeRB( this->vm()[
"sizeRB"].template as<int>() ),
105 muMin( this->vm()[
"muMin"].template as<double>() ),
106 muMax( this->vm()[
"muMax"].template as<double>() ),
107 mu( this->vm()[
"mu"].template as<double>() ),
108 Sampling( this->vm()[
"Sampling"].template as<int>() ),
109 SamplingCoarse( this->vm()[
"SamplingCoarse"].template as<int>() ),
110 Offline( this->vm()[
"Offline"].template as<int>() ),
111 ComputeError( this->vm()[
"ComputeError"].template as<int>() )
117 void run(
const double* X,
unsigned long P,
double* Y,
unsigned long N );
121 element_type blackbox( space_ptrtype Xh,
double param );
126 void ComputeSnapshot( space_ptrtype Xh, std::string filename );
128 void ChooseRBFunction( space_ptrtype Xh,
129 vector_of_element_type &VuBasis,
130 sparse_matrix_ptrtype
const &MassMatrix,
131 sparse_matrix_ptrtype
const &StiffMatrix );
133 void OrthogonalisationRBFunction( space_ptrtype Xh,
134 vector_of_element_type &Vu,
135 vector_of_element_type &M_Vu,
136 sparse_matrix_ptrtype
const &StiffMatrix,
137 sparse_matrix_ptrtype
const &MassMatrix );
139 void OrthogonalisationRBFunctionL2GrammSchmidt( space_ptrtype Xh,
140 vector_of_element_type &Vu,
141 sparse_matrix_ptrtype
const &MassMatrix );
143 double OrthogonalisationRBFunctionL2GrammSchmidt( space_ptrtype Xh,
144 vector_of_element_type &Vu,
146 sparse_matrix_ptrtype
const &MassMatrix );
149 Eigen::MatrixXd ConstructStiffMatrixSnapshot( space_ptrtype Xh
150 ,sparse_matrix_ptrtype
const &StiffMatrix );
152 Eigen::MatrixXd ConstructStiffMatrixRB( space_ptrtype Xh,
153 vector_of_element_type &Vu,
154 sparse_matrix_ptrtype
const &StiffMatrix );
156 Eigen::MatrixXd ConstructMassMatrixRB( space_ptrtype Xh,
157 vector_of_element_type &Vu,
158 sparse_matrix_ptrtype
const &MassMatrix );
160 void ConstructNIRB ( space_ptrtype Xh,
161 vector_of_element_type &M_VNirbBasis,
162 vector_of_element_type &VNirbBasis );
166 element_type BuildNirbSolutionWithoutPostProcess( space_ptrtype XhFine,
167 element_type uCoarseInterpolate,
168 vector_of_element_type
const &M_VNirbBasis,
169 vector_of_element_type
const &VNirbBasis,
170 Eigen::VectorXd & BetaiH );
171 element_type BuildCoarseInterpolation( space_ptrtype XhFine,
172 space_ptrtype XhCoarse,
175 element_type BuildNirbSolution( space_ptrtype XhFine,
176 space_ptrtype XhCoarse,
177 vector_of_element_type
const &M_VNirbBasis,
178 vector_of_element_type
const &VNirbBasis,
180 Eigen::VectorXd BetaiH );
183 Eigen::MatrixXd BuildBetaH( space_ptrtype XhFine,
184 space_ptrtype XhCoarse,
185 vector_of_element_type &M_VNirbBasis );
187 Eigen::MatrixXd BuildBetah( space_ptrtype XhFine,
188 vector_of_element_type &M_VNirbBasis );
190 Eigen::MatrixXd BuildPostProcessMatrix( space_ptrtype XhFine,
191 space_ptrtype XhCoarse,
192 vector_of_element_type &M_VNirbBasis );
194 element_type BuildNirbSolutionWithPostProcess( space_ptrtype XhFine,
195 space_ptrtype XhCoarse,
196 vector_of_element_type &M_VNirbBasis,
197 vector_of_element_type &VNirbBasis,
198 Eigen::VectorXd &BetaiuH,
199 element_type &uCoarseInterpolate );
203 gmsh_ptrtype createGeo(
double hsize,std::string MeshFileName );
208 backend_ptrtype M_backend;
211 double CoarseMeshSize;
216 int NbSnapshot,sizeRB;
217 double muMin,muMax,mu;
242 template<
int PolynomialOrder>
243 void NIRBTEST<PolynomialOrder>::run()
246 std::cout <<
"------------------------------------------------------------\n";
247 std::cout <<
"Execute NIRBTEST<" << PolynomialOrder <<
">\n";
248 std::vector<double> X( 12 );
250 X[1] = CoarseMeshSize;
251 X[2] = ReadingMeshes;
258 X[9] = SamplingCoarse;
260 X[11] = ComputeError;
261 std::vector<double> Y( 1 );
262 run( X.data(), X.size(), Y.data(), Y.size() );
266 template<
int PolynomialOrder>
267 void NIRBTEST<PolynomialOrder>::run(
const double* X,
unsigned long P,
double* Y,
unsigned long N )
270 std::cout<<
"Size reduced basis : " << sizeRB <<std::endl;
271 std::cout<<
"Polynomial order : P" << PolynomialOrder <<std::endl;
273 if ( !this->vm().count(
"nochdir" ) )
274 Environment::changeRepository( boost::format(
"%1%/P%2%/" )
275 % this->about().appName()
279 mesh_ptrtype meshExtraCoarse, meshCoarse, meshFine ;
283 meshCoarse = loadGMSHMesh( _mesh=
new mesh_type,
284 _filename=
"/home/chakir/Mesh/Mesh_H1.msh",
285 _update=MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
286 meshFine = loadGMSHMesh( _mesh=
new mesh_type,
287 _filename=
"/home/chakir/Mesh/Mesh_H2.msh",
288 _update=MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
289 std::cout <<
"Meshes read in file : done " <<std::endl;
294 std::cout <<
"Meshes build using 'createGMSHMesh' " <<std::endl;
295 std::cout <<
"Coarse mesh size : " << X[1] << std::endl;
296 meshCoarse = createGMSHMesh( _mesh=
new mesh_type,
297 _desc = createGeo( X[1],
"Mesh_Coarse" ),
303 std::cout <<
"Fine mesh size : " << X[0] << std::endl;
304 meshFine = createGMSHMesh( _mesh=
new mesh_type,
305 _desc = createGeo( X[0],
"Mesh_Fine" ),
312 std::cout <<
"WARNING : 'hfinesize > 'hcoarsesize' --> Fine mesh size set to hcoarsesize/2 : " << X[1]/2 << std::endl;
313 meshFine = createGMSHMesh( _mesh=
new mesh_type,
314 _desc = createGeo( X[1]/2,
"Mesh_Fine" ),
322 space_ptrtype XhFine = space_type::New( meshFine );
325 space_ptrtype XhCoarse = space_type::New( meshCoarse );
343 auto ui = XhFine->element();
345 vector_of_element_type VNirbBasis( sizeRB,ui );
346 vector_of_element_type MassMat_x_VuNirb( sizeRB,ui );
354 for (
int i = 0 ; i < sizeRB; i++ )
356 path = ( boost::format(
"./RB%1%NIRB_BasisFile_%2%" )%sizeRB%i ).str();
357 ui.load( _path=path );
359 if ( ui.l2Norm()== 0. )
361 std::cout <<
"WARNING Error : RB " << sizeRB <<
"NIRB_BasisFile_" << i <<
" is equal to zero => OFFline parameter set to 1 " <<std::endl;
367 path =( boost::format(
"./RB%1%D_Nirb_Basis%2%File" ) %sizeRB%i ).str() ;
368 ui.load( _path=path );
370 if ( ui.l2Norm()== 0. )
372 std::cout <<
"WARNING Error : RB" << sizeRB <<
"D_NIRB_Basis" << i <<
"File is equal to zero => Re-calculating M*Ui" <<std::endl;
373 auto Mui = M_backend->newVector( XhFine );
374 auto Ui = M_backend->newVector( XhFine );
376 auto uj = XhFine->element();
378 sparse_matrix_ptrtype MassMatrix = M_backend->newMatrix( _test=XhFine, _trial=XhFine );
379 form2( _test=XhFine , _trial=XhFine, _matrix= MassMatrix ) =
384 MassMatrix->multVector( Ui,Mui );
385 MassMat_x_VuNirb[i] = *Mui;
390 MassMat_x_VuNirb[i] = ui;
397 ConstructNIRB( XhFine,MassMat_x_VuNirb,VNirbBasis );
399 if ( SamplingCoarse )
401 std::cout <<
"Computation of Coarse snapshot " <<std::endl;
403 ComputeSnapshot( XhCoarse,
"_Coarse_" );
405 double Time_snapshot_Coarse = ti.elapsed();
406 std::cout <<
"Computation of the " << NbSnapshot <<
" snapshots : done -- ";
407 std::cout <<
"Time per snapshot: " << Time_snapshot_Coarse/NbSnapshot <<
" sec " <<std::endl;
420 std :: cout <<
"OFFLINE Procedure has been skipped" <<std::endl;
423 auto uCoarseInterpolate = XhFine->element();
424 auto uNirbCoarse = XhFine->element();
425 Eigen::VectorXd BetaiH( sizeRB );
426 double TimeCoarse,TimeCoarsePostProcess,TimeFine;
429 uCoarseInterpolate = BuildCoarseInterpolation( XhFine,XhCoarse,p );
430 TimeCoarse = ti.elapsed();
431 std :: cout <<
"Calculation of uCoarse and it's interpolation on XhFine :" << TimeCoarse <<
" sec" <<std::endl;
432 uNirbCoarse = BuildNirbSolutionWithoutPostProcess( XhFine,uCoarseInterpolate,MassMat_x_VuNirb,VNirbBasis,BetaiH );
434 TimeCoarse = ti.elapsed();
435 std :: cout <<
"Construction of NIRB solution (uNirbCoarse) - Fine/Coarse Grid (saved in nirb2GridCoarse):" <<std::endl;
436 std::cout <<
"Time to build solution " << TimeCoarse <<
" sec" <<std::endl;
438 export_ptrtype exporter2GridCoarse( export_type::New( this->vm(),
"nirb2GridCoarse" ) );
439 exporter2GridCoarse->step( 0 )->setMesh( meshFine );
440 exporter2GridCoarse->step( 0 )->add(
"uNirbCoarse",uNirbCoarse );
441 exporter2GridCoarse->save();
443 auto uNirbCoarsePostProcess = XhFine->element();
447 std::cout <<
"Computation of Coarse snapshot " <<std::endl;
449 ComputeSnapshot( XhCoarse,
"_Coarse_" );
451 double Time_snapshot_Coarse = ti.elapsed();
452 std::cout <<
"Computation of the " << NbSnapshot <<
" snapshots : done -- ";
453 std::cout <<
"Time per snapshot: " << Time_snapshot_Coarse/NbSnapshot <<
" sec " <<std::endl;
455 uNirbCoarsePostProcess = BuildNirbSolutionWithPostProcess( XhFine, XhCoarse, MassMat_x_VuNirb,
456 VNirbBasis,BetaiH,uCoarseInterpolate );
458 TimeCoarsePostProcess = ti.elapsed();
459 std :: cout <<
"Post-processing of NIRB solution (uNirbCoarsePostProcess) - Fine/Coarse Grid (saved in nirb2GridCoarse):" <<std::endl;
460 std::cout <<
"Time " << TimeCoarsePostProcess <<
" sec" <<std::endl;
462 export_ptrtype exporter2GridCoarsePostProcess( export_type::New( this->vm(),
"nirb2GridCoarsePostProcess" ) );
463 exporter2GridCoarsePostProcess->step( 0 )->setMesh( meshFine );
464 exporter2GridCoarsePostProcess->step( 0 )->add(
"uNirbCoarsePostProcess",uNirbCoarsePostProcess );
465 exporter2GridCoarsePostProcess->save();
469 std :: cout <<
"Error calculation " <<std::endl;
470 Eigen::VectorXd Betaih( sizeRB );
472 auto u1Grid = XhFine->element();
474 u1Grid = BuildNirbSolution( XhFine,XhFine,MassMat_x_VuNirb,VNirbBasis,p,Betaih );
475 TimeFine = ti.elapsed();
476 std::cout <<
"Construction of uNirbFine - Fine/Fine Grid (saved in nirb1Grid): done "<<std::endl;
477 std::cout <<
"Time to build " << TimeFine <<
" sec" <<std::endl;
481 mesh_ptrtype meshRef;
485 meshRef = loadGMSHMesh( _mesh=
new mesh_type,
486 _filename=
"/home/chakir/Mesh/Mesh_H3.msh",
487 _update=MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
494 meshRef = createGMSHMesh( _mesh=
new mesh_type,
495 _desc = createGeo( X[0]/2,
"Mesh_Fine" ),
501 meshRef = createGMSHMesh( _mesh=
new mesh_type,
502 _desc = createGeo( X[1]/4,
"Mesh_Fine" ),
509 space_ptrtype XhRef = space_type::New( meshRef );
517 auto uRef = XhRef->element();
519 uRef = blackbox( XhRef, p );
520 std :: cout <<
"Computation of FE solution (uRef) - Ref Grid (not saved): done" <<std::endl;
525 _expr=( idv( uRef )*idv( uRef ) ) ).evaluate()( 0,0 ) );
526 double SemiH1NormUref = std :: sqrt(
integrate( _range=
elements( XhRef->mesh() ),
527 _expr=( gradv( uRef )*trans( gradv( uRef ) ) ) ).evaluate()( 0,0 ) );
528 double H1NormUref = std :: sqrt( SemiH1NormUref *SemiH1NormUref + L2NormUref*L2NormUref );
548 double ErrH1uNirbCoarse,ErrH1uNirbCoarsePostProcess,ErrH1u1Grid;
554 _expr=( ( gradv( uRef )-gradv( uNirbCoarsePostProcess ) )*trans( gradv( uRef )-gradv( uNirbCoarsePostProcess ) )
555 + ( idv( uRef )-idv( uNirbCoarsePostProcess ) )*( idv( uRef )-idv( uNirbCoarsePostProcess ) ) ) ).evaluate()( 0,0 );
556 ErrH1uNirbCoarsePostProcess = sqrt( ErrH1uNirbCoarsePostProcess )/H1NormUref;
561 _expr=( ( gradv( uRef )-gradv( uNirbCoarse ) )*trans( gradv( uRef )-gradv( uNirbCoarse ) )
562 + ( idv( uRef )-idv( uNirbCoarse ) )*( idv( uRef )-idv( uNirbCoarse ) ) ) ).evaluate()( 0,0 );
563 ErrH1uNirbCoarse = sqrt( ErrH1uNirbCoarse )/H1NormUref;
567 _expr=( ( gradv( uRef )-gradv( u1Grid ) )*trans( gradv( uRef )-gradv( u1Grid ) )
568 + ( idv( uRef )-idv( u1Grid ) )*( idv( uRef )-idv( u1Grid ) ) ) ).evaluate()( 0,0 );
569 ErrH1u1Grid = sqrt( ErrH1u1Grid )/H1NormUref;
571 std :: cout <<
"H1-norm NIRB Error " <<std::endl;
572 std :: cout <<
"||u_ref - u_NirbCoarse ||_{H1} = " << ErrH1uNirbCoarse <<std::endl;
573 std :: cout <<
"||u_ref - u_NirbCoarsePostProcess ||_{H1} = " << ErrH1uNirbCoarsePostProcess <<std::endl;
574 std :: cout <<
"||u_ref - u_NirbFine ||_{H1} = " << ErrH1u1Grid<<std::endl;
575 std :: cout << std::endl;
581 export_ptrtype exporter1Grid( export_type::New( this->vm(),
"nirbFine" ) );
582 exporter1Grid->step( 0 )->setMesh( meshFine );
583 exporter1Grid->step( 0 )->add(
"uNirbFine",u1Grid );
584 exporter1Grid->save();
588 auto uFine = XhFine->element();
589 auto uCoarse = XhCoarse->element();
591 uCoarse = blackbox( XhCoarse, p );
592 std :: cout <<
"Computation of FE solution (uCoarse) - Coarse Grid (not saved):"<<std::endl;
594 uFine = blackbox( XhFine, p );
595 std :: cout <<
"Computation of FE solution (uFine) - Fine Grid (saved in EFFine): done" <<std::endl;
637 export_ptrtype exporterFine( export_type::New( this->vm(),
"EF_Fine" ) );
638 export_ptrtype exporterCoarse( export_type::New( this->vm(),
"EF_Coarse" ) );
639 exporterFine->step( 0 )->setMesh( meshFine );
640 exporterCoarse->step( 0 )->setMesh( meshCoarse );
642 exporterFine->step( 0 )->add(
"uFine", uFine );
643 exporterCoarse->step( 0 )->add(
"uCoarse", uCoarse );
645 exporterFine->save();
646 exporterCoarse->save();
662 template<
int PolynomialOrder>
663 gmsh_ptrtype NIRBTEST<PolynomialOrder>::createGeo(
double hsize, std::string MeshFileName )
665 double hCornersize1 = hsize/15.;
666 double hCornersize2 = hsize/30.;
667 std::ostringstream ostr;
668 ostr <<
"Mesh.MshFileVersion = 2.2;" <<
"\n";
669 ostr <<
"Mesh.CharacteristicLengthExtendFromBoundary=1;" <<std::endl;
670 ostr <<
"Mesh.CharacteristicLengthFromPoints=1;" <<std::endl;
671 ostr <<
"Mesh.ElementOrder=1;" <<std::endl;
672 ostr <<
"Mesh.SecondOrderIncomplete = 0;" <<std::endl;
673 ostr <<
"Mesh.Algorithm = 6;" <<std::endl;
674 ostr <<
"Mesh.OptimizeNetgen=1;" <<std::endl;
675 ostr <<
"// partitioning data" <<std::endl;
676 ostr <<
"Mesh.Partitioner=1;" <<std::endl;
677 ostr <<
"Mesh.NbPartitions=1;" <<std::endl;
678 ostr <<
"Mesh.MshFilePartitioned=0;" <<std::endl;
679 ostr <<
"Point(1) = {0,0,0.0," << hsize <<
"};" <<std::endl;
680 ostr <<
"Point(2) = {1,0,0.0," << hCornersize2 <<
"};" <<std::endl;
681 ostr <<
"Point(3) = {1,1,0.0," << hCornersize1 <<
"};"<<std::endl;
682 ostr <<
"Point(4) = {0,1,0.0," << hCornersize2 <<
"};"<<std::endl;
683 ostr <<
"Line(1) = {4,1};" <<std::endl;
684 ostr <<
"Line(2) = {1,2};" <<std::endl;
685 ostr <<
"Line(3) = {2,3};" <<std::endl;
686 ostr <<
"Line(4) = {3,4};" <<std::endl;
687 ostr <<
"Line Loop(5) = {1,2,3,4};" <<std::endl;
688 ostr <<
"Plane Surface(6) = {5};" <<std::endl;
689 ostr <<
"Physical Line(\"Dirichlet1\") = {1,2};" <<std::endl;
690 ostr <<
"Physical Line(\"Dirichlet2\") = {3};" <<std::endl;
691 ostr <<
"Physical Line(\"Dirichlet3\") = {4};" <<std::endl;
692 ostr <<
"Physical Surface(\"Mat1\") = {6};" <<std::endl;
693 std::ostringstream nameStr;
694 nameStr.precision( 3 );
695 nameStr << MeshFileName;
697 gmsh_ptrtype gmshp(
new Gmsh );
698 gmshp->setPrefix( nameStr.str() );
699 gmshp->setDescription( ostr.str() );
707 template<
int PolynomialOrder>
708 void NIRBTEST<PolynomialOrder>::ComputeSnapshot( space_ptrtype Xh,std::string filename )
711 int XhNdof = Xh->nLocalDof();
712 double pi = 3.14159265;
714 auto ui = Xh->element();
716 for (
int i =0; i<NbSnapshot; i++ )
718 double theta = ( i+1 )*( muMax-muMin )/( 113. );
719 ui = blackbox( Xh,theta );
722 std::cerr <<
"In ComputeSnapshot : ERROR IN USING BLACKBOX -> ui = 0" << endl;
725 std::string path =
"./Sol" + filename + ( boost::format(
"%1%" ) %i ).str() ;
726 ui.save( _path=path );
736 template<
int PolynomialOrder>
737 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: ConstructStiffMatrixSnapshot( space_ptrtype Xh,sparse_matrix_ptrtype
const & StiffMatrix )
741 auto ui = Xh->element();
742 auto uj = Xh->element();
748 Eigen::MatrixXd S ( NbSnapshot,NbSnapshot );
749 cout <<
" Dans ConstructStiffMatrixSnapshot " << endl;
753 for (
int i = 0; i< NbSnapshot; i++ )
756 std::string path = ( boost::format(
"./Sol_%1%" ) %i ).str() ;
760 ui.load( _path=path);
762 if ( ui.l2Norm()==0. )
764 std::cerr <<
"In 'ConstructStiffMatrixSnapshot': ERROR IN LOADING FILE " << path <<std::endl;
765 auto t = Xh->element();
767 cout <<
"Reloading of ui done -> L2norm(ui)= " << t.l2Norm() << endl;
771 for (
int j=0; j< i; j++ )
774 path = ( boost::format(
"./Sol_%1%" ) %j ).str() ;
775 uj.load( _path=path );
777 if ( uj.l2Norm()==0. )
779 std::cerr <<
"In 'ConstructStiffMatrixSnapshot': ERROR IN LOADING FILE " << path <<std::endl;
783 double Sij = StiffMatrix->energy( ui,uj );
788 double Sii = StiffMatrix->energy( ui,ui );
792 std::string filename = ( boost::format(
"StiffMatrixP%1%" )%PolynomialOrder ).str();
793 std::ofstream fileMat( filename.c_str() );
794 fileMat << NbSnapshot <<std::endl;
796 for (
int i = 0; i<NbSnapshot; i++ )
798 for (
int j =0; j<NbSnapshot; j++ )
800 fileMat << S( i,j ) <<std::endl;
812 template<
int PolynomialOrder>
813 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: ConstructStiffMatrixRB( space_ptrtype Xh,vector_of_element_type &Vu,sparse_matrix_ptrtype
const & StiffMatrix )
816 Eigen::MatrixXd S ( sizeRB,sizeRB );
818 for (
int i = 0; i< sizeRB; i++ )
820 for (
int j=0; j< i; j++ )
822 double Sij = StiffMatrix->energy( Vu[i],Vu[j] );
827 double Sii = StiffMatrix->energy( Vu[i],Vu[i] );
835 template<
int PolynomialOrder>
836 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: ConstructMassMatrixRB( space_ptrtype Xh,vector_of_element_type &Vu,sparse_matrix_ptrtype
const &MassMatrix )
839 Eigen::MatrixXd S ( sizeRB,sizeRB );
841 for (
int i = 0; i< sizeRB; i++ )
843 for (
int j=0; j< i; j++ )
845 double Sij = MassMatrix->energy( Vu[i],Vu[j] );
850 double Sii = MassMatrix->energy( Vu[i],Vu[i] );
860 template<
int PolynomialOrder>
861 void NIRBTEST<PolynomialOrder> ::ChooseRBFunction( space_ptrtype Xh,vector_of_element_type &VuBasis,sparse_matrix_ptrtype
const &MassMatrix,sparse_matrix_ptrtype
const &StiffMatrix )
865 Eigen::MatrixXd S ( NbSnapshot,NbSnapshot );
869 std::string str = ( boost::format(
"StiffMatrixP%1%" )%PolynomialOrder ).str();
870 std::ifstream fileMat( str.c_str() );
874 std::cerr <<
"Error in reading file " << ( boost::format(
"StiffMatrixP%1%" )%PolynomialOrder ).str() <<
"does not exist " << std::endl;
875 std::cout <<
"Need to re-build the Stiffness Matrix -> SAMPLING SET TO 1 " << std::endl;
884 if ( Itemp != NbSnapshot )
886 std :: cerr <<
"Error in reading file" << ( boost::format(
"StiffMatrixP%1%" )%PolynomialOrder ).str() <<
" -- " << Itemp <<
" != NbSnapshot(" << NbSnapshot <<
")" << std::endl;
887 std::cout <<
"Need to re-build the Stiffness Matrix -> SAMPLING SET TO 1 " << std::endl;
893 for (
int i=0; i<NbSnapshot; i++ )
895 for (
int j =0; j<NbSnapshot; j++ )
909 S = ConstructStiffMatrixSnapshot( Xh,StiffMatrix );
913 Eigen::EigenSolver <Eigen::MatrixXd> eigen_solver( S );
918 std :: cout <<
"Computation of the eigenvalues of the stiffness matrix S (NbSnapshot,NbSnapshot) : done" <<std::endl;
927 std::string path = ( boost::format(
"./IndBR%1%" ) %sizeRB ).str() ;
928 std::ofstream find( path.c_str() );
932 std :: cerr <<
" 'ChooseRBFunction routine' - ERROR IN OPENING FILE:" << path <<std::endl;
936 std::vector<int> Uind( NbSnapshot );
938 for (
int i=0; i<NbSnapshot; i++ )
944 auto ui = Xh->element();
950 std::vector<double>Vi( NbSnapshot );
952 for (
int i=0; i<sizeRB; i++ )
954 for (
int ri = 0; ri < NbSnapshot; ri++ )
960 Vi[ri] = std::abs( real( eigen_solver.eigenvectors().col( i )[ri] ) );
971 int IndMax = std::distance( Vi.begin(),boost::max_element( Vi ) );
975 path = ( boost::format(
"./Sol_%1%" ) %IndMax ).str();
976 ui.load( _path=path );
979 if ( ui.l2Norm()==0. )
981 std::cerr <<
"'ChooseRBFunction routine': ERROR IN LOADING FILE " << path <<std::endl;
989 double normL2 = OrthogonalisationRBFunctionL2GrammSchmidt( Xh,VuBasis,i,MassMatrix );
996 find << IndMax <<std::endl;
1009 find << IndMax <<std::endl;
1021 for (
int i=0; i<sizeRB; i++ )
1023 L2normVui = std::sqrt( MassMatrix->energy( VuBasis[i],VuBasis[i] ) );
1024 VuBasis[i].scale( 1./L2normVui );
1025 VuBasis[i].save( _path=( boost::format(
"./Sol_OR%1%" ) %i ).str() );
1034 template<
int PolynomialOrder>
1035 double NIRBTEST <PolynomialOrder> :: OrthogonalisationRBFunctionL2GrammSchmidt( space_ptrtype Xh, vector_of_element_type &Vu,
int n,sparse_matrix_ptrtype
const &MassMatrix )
1038 if ( MassMatrix->linftyNorm()==0. )
1040 std::cout <<
"ERROR in OrthogonalisationRBFunctionL2GrammSchmidt : MassMatrix is null " <<std::endl;
1044 double Dtemp1,Dtemp2,Dtemp3,Dtemp4;
1047 for (
int k=0; k<n; k++ )
1050 Dtemp1 = MassMatrix->energy( Vu[n],Vu[k] );
1051 Dtemp2 = MassMatrix->energy( Vu[k],Vu[k] );
1052 Dtemp3 = -Dtemp1/Dtemp2;
1053 Vu[n].add( Dtemp3,Vu[k] );
1057 Dtemp4 = MassMatrix->energy( Vu[n],Vu[n] );
1065 template<
int PolynomialOrder>
1066 void NIRBTEST<PolynomialOrder> :: OrthogonalisationRBFunctionL2GrammSchmidt
1067 ( space_ptrtype Xh,vector_of_element_type &Vu, sparse_matrix_ptrtype
const &MassMatrix )
1070 double Dtemp,Dtemp1,Dtemp2,Dtemp3, L2norm;
1072 for (
int i=1; i<sizeRB; i++ )
1074 for (
int k=0; k<i; k++ )
1076 Dtemp1 = MassMatrix->energy( Vu[i],Vu[k] );
1077 Dtemp2 = MassMatrix->energy( Vu[k],Vu[k] );
1078 Dtemp3 = - Dtemp1/Dtemp2;
1079 Vu[i].add( Dtemp3,Vu[k] );
1084 for (
int i=0; i<sizeRB; i++ )
1086 L2norm = std::sqrt( MassMatrix->energy( Vu[i],Vu[i] ) );
1087 Vu[i].scale( 1./L2norm ) ;
1094 template<
int PolynomialOrder>
1095 void NIRBTEST<PolynomialOrder> ::OrthogonalisationRBFunction( space_ptrtype Xh,vector_of_element_type &Vu,vector_of_element_type &M_Vu,sparse_matrix_ptrtype
const &StiffMatrix,sparse_matrix_ptrtype
const &MassMatrix )
1100 auto ui = Xh->element();
1111 Eigen :: MatrixXd A ( sizeRB,sizeRB );
1112 Eigen :: MatrixXd B ( sizeRB,sizeRB );
1113 Eigen :: MatrixXd C ( sizeRB,sizeRB );
1115 A = ConstructStiffMatrixRB( Xh,Vu,StiffMatrix );
1116 B = ConstructMassMatrixRB( Xh,Vu,MassMatrix );
1120 Eigen:: EigenSolver<Eigen::MatrixXd> es;
1123 vector_of_element_type VuNirb( sizeRB,ui );
1126 for (
int i =0; i<sizeRB; i++ )
1128 for (
int k=0; k<sizeRB; k++ )
1130 VuNirb[i].add( real( es.eigenvectors().col( i )[k] ),Vu[k] );
1135 OrthogonalisationRBFunctionL2GrammSchmidt( Xh,VuNirb,MassMatrix );
1137 auto Mui = M_backend->newVector( Xh );
1138 auto Ui = M_backend->newVector( Xh );
1141 for (
int i=0; i<sizeRB; i++ )
1143 path = ( boost::format(
"./RB%1%NIRB_BasisFile_%2%" ) %sizeRB%i ).str() ;
1144 VuNirb[i].save( _path=path );
1148 MassMatrix->multVector( Ui,Mui );
1151 if ( M_Vu[i].l2Norm()>1e-20 )
1153 path = ( boost::format(
"./RB%1%D_Nirb_Basis%2%File" )%sizeRB%i ).str();
1154 M_Vu[i].save( _path=path );
1159 std::cerr <<
"ERROR: in computation of D*NIRB Basis function (" << i <<
") : Result equal to zero " <<std::endl;
1169 template<
int PolynomialOrder>
1170 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder> :: BuildNirbSolutionWithoutPostProcess ( space_ptrtype XhFine, element_type uCoarseInterpolate,vector_of_element_type
const &M_VNirbBasis, vector_of_element_type
const &VNirbBasis, Eigen::VectorXd & BetaiH )
1175 auto uNirb = XhFine->element();
1177 auto UcoarseInterpol = M_backend->newVector( XhFine );
1178 *UcoarseInterpol = uCoarseInterpolate;
1184 auto Mui = M_backend->newVector( XhFine );
1186 for (
int i =0; i<sizeRB; i++ )
1188 *Mui=M_VNirbBasis[i];
1190 uNirb.add( BetaiH[i], VNirbBasis[i] );
1201 template<
int PolynomialOrder>
1202 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder> ::BuildCoarseInterpolation( space_ptrtype XhFine,space_ptrtype XhCoarse,
double param )
1206 auto uCoarse = XhCoarse->element();
1207 auto uCoarseInterpolate = XhFine->element();
1208 auto ui = XhFine->element();
1209 uCoarse = blackbox( XhCoarse,param );
1219 interpolate ( XhFine,uCoarse,uCoarseInterpolate );
1220 return uCoarseInterpolate;
1228 template<
int PolynomialOrder>
1229 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder> :: BuildNirbSolution( space_ptrtype XhFine,space_ptrtype XhCoarse, vector_of_element_type
const &M_VNirbBasis, vector_of_element_type
const &VNirbBasis,
double param,Eigen::VectorXd BetaiH )
1233 auto uCoarseInterpolate = XhFine->element();
1234 auto uNirb = XhFine->element();
1235 uCoarseInterpolate = BuildCoarseInterpolation( XhFine,XhCoarse,param );
1236 uNirb = BuildNirbSolutionWithoutPostProcess( XhFine,uCoarseInterpolate,M_VNirbBasis,VNirbBasis,BetaiH );
1242 template<
int PolynomialOrder>
1243 void NIRBTEST<PolynomialOrder> :: ConstructNIRB( space_ptrtype Xh, vector_of_element_type &M_VNirbBasis, vector_of_element_type &VNirbBasis )
1246 std :: cout <<
"OFFLINE PROCEDURE : Construction of the 'non intruisive' reduced basis (nirb) " <<std::endl;
1248 double Time_snapshot = 0.;
1252 std::cout <<
"Sampling Procedure :" <<std::endl;
1253 ComputeSnapshot( Xh,
"_" );
1254 Time_snapshot = ti.elapsed();
1255 std::cout <<
"Computation of the " << NbSnapshot <<
" snapshots : done -- ";
1256 std::cout <<
"Time per snapshot: " << Time_snapshot/NbSnapshot <<
" sec " <<std::endl;
1262 auto ui = Xh->element();
1263 auto uj = Xh->element();
1265 sparse_matrix_ptrtype MassMatrix = M_backend->newMatrix( _test=Xh, _trial=Xh );
1266 form2( _test=Xh , _trial=Xh, _matrix= MassMatrix ) =
1269 sparse_matrix_ptrtype StiffMatrix = M_backend->newMatrix( _test=Xh, _trial=Xh );
1270 form2( _test=Xh , _trial=Xh, _matrix= StiffMatrix ) =
1271 integrate( _range=
elements( Xh->mesh() ), _expr=gradt( ui )*trans( grad( uj ) ) );
1272 double Time_Construct_Elementary_Matrix = ti.elapsed();
1273 std::cout <<
"Time to build elementary F.E matrix = " << Time_Construct_Elementary_Matrix <<
" sec " <<std::endl;
1276 ChooseRBFunction( Xh,VNirbBasis,MassMatrix,StiffMatrix );
1277 double TimeChooseRB = ti.elapsed();
1278 std::cout <<
"Choice of " << sizeRB <<
" reduced basis functions : done -- ";
1279 std::cout <<
"Time: " << TimeChooseRB <<
" sec "<<std::endl;
1284 OrthogonalisationRBFunction( Xh,VNirbBasis,M_VNirbBasis,StiffMatrix,MassMatrix );
1286 double Time_buildNIRBbasis = ti.elapsed();
1287 std::cout <<
"H1 and L2 orthogonalisation of the reduced basis functions : done -- ";
1288 std::cout <<
"Time: " << Time_buildNIRBbasis <<
" sec " <<std::endl;
1289 double TotalTime = Time_buildNIRBbasis +TimeChooseRB + Time_Construct_Elementary_Matrix + Time_snapshot;
1290 std::cout <<
"Total time for the OFFLINE procedure: " << TotalTime <<
" sec "<<std::endl;
1300 template<
int PolynomialOrder>
1301 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder>::blackbox( space_ptrtype Xh,
double param )
1304 auto u = Xh->element();
1305 auto v = Xh->element();
1307 value_type pi = M_PI;
1308 double penaltyTerm = 30.;
1311 auto velocity = vec( cst( cos( param ) ),cst( sin( param ) ) );
1317 auto g = ( Px()*Px()*Py()*Py() );
1318 auto F = M_backend->newVector( Xh );
1324 auto D = M_backend->newMatrix( _test=Xh, _trial=Xh );
1325 form2( _test=Xh, _trial=Xh, _matrix=D ) =
1326 integrate( _range=
elements( Xh->mesh() ), _expr=0.01*gradt( u )*trans( grad( v ) )+ ( gradt( u )*velocity )*
id( v ) );
1334 form2( _test=Xh, _trial=Xh, _matrix=D ) += on(
boundaryfaces( Xh->mesh() ), u, F,g );
1344 template<
int PolynomialOrder>
1345 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: BuildBetaH( space_ptrtype XhFine,
1346 space_ptrtype XhCoarse,vector_of_element_type &M_VNirbBasis )
1354 auto uCoarse = XhCoarse->element();
1355 auto uCoarseInterpolation = XhFine->element();
1357 std::vector<int> IndTab( sizeRB );
1358 Eigen::MatrixXd BetaH( sizeRB,sizeRB );
1359 std::string path = ( boost::format(
"./IndBR%1%" ) %sizeRB ).str() ;
1360 std::ifstream find( path.c_str() );
1364 std::cerr <<
"ERROR : Problem in opening 'IndBR'" << sizeRB <<
" file " <<std::endl;
1370 for (
int i =0 ; i < sizeRB; i++ )
1377 path = ( boost::format(
"./AlphaH_%1%" ) %sizeRB ).str();
1378 std::ofstream fB( path.c_str() );
1379 fB << sizeRB <<
" " << sizeRB <<std::endl;
1381 for (
int i =0; i<sizeRB; i++ )
1384 uCoarseInterpolation.zero();
1385 path = ( boost::format(
"./Sol_Coarse_%1%" ) %IndTab[i] ).str();
1386 uCoarse.load( _path=path );
1388 if ( uCoarse.l2Norm()==0. )
1390 std::cerr <<
"'BuildBetaH routine': ERROR IN LOADING FILE " << path <<std::endl;
1391 std::cerr <<
"Coarse sampling has to be done " << std::endl;
1392 ComputeSnapshot( XhCoarse,
"_Coarse_" );
1398 interpolate( XhFine,uCoarse,uCoarseInterpolation );
1399 auto UcoarseInterpol_i = M_backend->newVector( XhFine );
1400 *UcoarseInterpol_i = uCoarseInterpolation;
1401 auto Muj = M_backend->newVector( XhFine );
1403 for (
int j =0; j<sizeRB; j++ )
1405 *Muj=M_VNirbBasis[j];
1407 fB << BetaH( j,i ) <<std::endl;
1418 template<
int PolynomialOrder>
1419 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: BuildBetah( space_ptrtype XhFine,
1420 vector_of_element_type &M_VNirbBasis )
1423 auto uFine = XhFine->element();
1425 std::vector<int> IndTab( sizeRB );
1426 Eigen::MatrixXd BetaH( sizeRB,sizeRB );
1427 std::string path = ( boost::format(
"./IndBR%1%" ) %sizeRB ).str() ;
1428 std::ifstream find( path.c_str() );
1432 std::cerr <<
"ERROR : Problem in opening 'IndBR'" << sizeRB <<
" file " <<std::endl;
1436 for (
int i =0 ; i < sizeRB; i++ )
1443 path = ( boost::format(
"./Betah_%1%" ) %sizeRB ).str();
1444 std::ofstream fB( path.c_str() );
1445 fB << sizeRB <<
" " << sizeRB <<std::endl;
1447 for (
int i =0; i<sizeRB; i++ )
1450 path = ( boost::format(
"./Sol_%1%" ) %IndTab[i] ).str();
1451 uFine.load( _path=path );
1453 if ( uFine.l2Norm()==0. )
1455 std::cerr <<
"'ChooseRBFunction routine': ERROR IN LOADING FILE " << path <<std::endl;
1459 auto Ufine_i = M_backend->newVector( XhFine );
1461 auto Muj = M_backend->newVector( XhFine );
1463 for (
int j =0; j<sizeRB; j++ )
1465 *Muj=M_VNirbBasis[j];
1467 fB << BetaH( j,i ) <<std::endl;
1479 template<
int PolynomialOrder>
1480 Eigen::MatrixXd NIRBTEST<PolynomialOrder> ::BuildPostProcessMatrix( space_ptrtype XhFine, space_ptrtype XhCoarse,vector_of_element_type &M_VNirbBasis )
1483 Eigen::MatrixXd BetaH( sizeRB,sizeRB );
1484 Eigen::MatrixXd Betah( sizeRB,sizeRB );
1485 Eigen::MatrixXd T( sizeRB,sizeRB );
1487 int FileOkA = 1,FileOkB = 1 ;
1491 std::string path = ( boost::format(
"./AlphaH_%1%" ) %sizeRB ).str();
1492 std::ifstream fA ( path.c_str() );
1497 BetaH = BuildBetaH( XhFine,XhCoarse,M_VNirbBasis );
1500 path = ( boost::format(
"./betah_%1%" ) %sizeRB ).str();
1501 std::ifstream fB ( path.c_str() );
1506 Betah = BuildBetah( XhFine,M_VNirbBasis );
1510 double dTempA,dTempB;
1512 if ( FileOkA && FileOkB )
1514 fA >> Itemp1 >> Itemp2;
1515 fB >> Itemp1 >> Itemp2;
1517 for (
int i = 0; i< sizeRB; i++ )
1519 for (
int j=0; j<sizeRB; j++ )
1523 BetaH( j,i ) = dTempA;
1524 Betah( j,i ) = dTempB;
1533 fA >> Itemp1 >> Itemp2;
1535 for (
int i = 0; i< sizeRB; i++ )
1537 for (
int j=0; j<sizeRB; j++ )
1540 BetaH( j,i ) = dTempA;
1547 fB >> Itemp1 >> Itemp2;
1549 for (
int i = 0; i< sizeRB; i++ )
1551 for (
int j=0; j<sizeRB; j++ )
1554 Betah( j,i ) = dTempB;
1567 BetaH = BuildBetaH( XhFine,XhCoarse,M_VNirbBasis );
1568 Betah = BuildBetah( XhFine,M_VNirbBasis );
1571 T = Betah* BetaH.inverse() ;
1582 template<
int PolynomialOrder>
1583 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder> ::
1584 BuildNirbSolutionWithPostProcess( space_ptrtype XhFine,
1585 space_ptrtype XhCoarse,
1586 vector_of_element_type &M_VNirbBasis,
1587 vector_of_element_type &VNirbBasis,
1588 Eigen::VectorXd &BetaiuH,
1589 element_type &uCoarseInterpolate )
1592 Eigen::MatrixXd T ( sizeRB,sizeRB );
1593 T = BuildPostProcessMatrix( XhFine,XhCoarse,M_VNirbBasis );
1598 auto uNirb = XhFine->element();
1601 for (
int i =0; i<sizeRB; i++ )
1605 for (
int k=0; k<sizeRB; k++ )
1607 AlphaiH += T( i,k )*BetaiuH[k];
1610 uNirb.add( AlphaiH, VNirbBasis[i] );