GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemReducedSteklov.hxx
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 111 304 36.5%
Branches: 105 659 15.9%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors:
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "FiniteElement/elementMatrix.hpp"
21 #include "Core/felisceTransient.hpp"
22 #include "Geometry/curvatures.hpp"
23 #include "Solver/steklovBanner.hpp"
24
25 namespace felisce {
26
27 // ======================================================================
28 // the public methods
29 // ======================================================================
30 /* \brief Constructor
31 */
32 template <class volumeProblem>
33 4 LinearProblemReducedSteklov<volumeProblem>::LinearProblemReducedSteklov():
34 volumeProblem(),
35 4 m_reducedSteklov(nullptr),
36
4/8
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 17 taken 4 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 4 times.
✗ Branch 21 not taken.
 4 m_chronoRS(nullptr)
37 {
38 4 m_imgType=dirichlet; // Default value is Dirichlet, but it can be changed in derived classes.
39 4 }
40
41 /*! \brief Destructor
42 *
43 * We destroy the reducedSteklov class
44 */
45 template <class volumeProblem>
46 8 LinearProblemReducedSteklov<volumeProblem>::~LinearProblemReducedSteklov(){
47 8 if ( m_reducedSteklov )
48 ✗ delete m_reducedSteklov;
49 }
50
51 /* \brief Function to assemble the system of the volume problem
52 * Matrix and Boundary Conditions
53 */
54 template <class volumeProblem>
55 void
56 656 LinearProblemReducedSteklov<volumeProblem>::assembleVolumeSystem( FlagMatrixRHS flag ) {
57 656 this->assembleMatrixRHS(MpiInfo::rankProc(), flag );
58 656 this->derivedProblemAssemble(flag);
59
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 328 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 656 if ( m_imgType == neumann && !( flag == FlagMatrixRHS::only_rhs ) ) {
60 ✗ this->createAndCopyMatrixRHSWithoutBC();
61 }
62 656 this->finalizeEssBCTransient();
63 656 this->applyBC( FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), flag, flag );
64 }
65 /*! \brief Computation/Loading of the real (Full rank) Steklov operator
66 *
67 * this function compute the Steklov-Poincare operator
68 * the idea is to solve N times the volume problem
69 * where N is the dimension of the interface
70 * each time the data is computed via dirac delta in a different dof
71 */
72 template <class volumeProblem>
73 void
74 4 LinearProblemReducedSteklov<volumeProblem>::assembleFullRankSteklov() {
75 4 this->initChronoRS();
76 //! We allocate the (dense) matrix for the Steklov operator
77
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
78 4 m_fullSteklov.createDense(this->m_dofBD[0/*iBD*/].comm(),
79 4 /*local row size*/ this->m_dofBD[0/*iBD*/].numLocalDofInterface(),
80 4 /*local col size*/ this->m_dofBD[0/*iBD*/].numLocalDofInterface(),
81 4 /*global row size*/ this->m_dofBD[0/*iBD*/].numGlobalDofInterface(),
82 4 /*global col size*/ this->m_dofBD[0/*iBD*/].numGlobalDofInterface(),
83 /*where to allocate the matrix*/FELISCE_PETSC_NULLPTR);
84 4 m_fullSteklov.setFromOptions();
85 }
86
87 //! The operator is computed via the solution of the volume problem.
88 //! this operator is obviously skipped when the operator is loaded from file.
89
1/3
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 switch ( FelisceParam::instance().optionForSteklov ) {
90 ✗ case LOAD_IT:
91 ✗ break;
92 4 case COMPUTE_IT:
93 case COMPUTE_AND_SAVE_IT:
94 {
95 // we raise this flag to use correctly the boundary conditions
96 // function
97
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 this->useSteklovDataBegin();
98
99 //! We use a banner object to handle the comunications on the screen of this function
100 //! and to improve readibility
101
2/4
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
 4 steklovBanner banner(this->m_dofBD[0/*iBD*/].numLocalDofInterface(), MpiInfo::rankProc(), this->m_dofBD[0/*iBD*/].comm());
102
103
104 // We switch off the verbosity of the KSP
105
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 FelisceParam::linearSolverVerboseOFF();
106
107 // We assemble the matrix of the linearProblem ( the volume matrix )
108 // including the informations needed for the BC
109 // the flag is std::set to only_matrix since we want to build just the matrix
110
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 this->clearMatrixRHS( FlagMatrixRHS::only_matrix );
111
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 this->assembleVolumeSystem( FlagMatrixRHS::only_matrix );
112
113
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 banner.initComputation();
114 std::size_t idDofGlobalApplication, idDofGlobalVolumePetsc;
115
3/4
✓ Branch 1 taken 324 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 324 times.
✓ Branch 6 taken 4 times.
 328 for ( idDofGlobalApplication=0; idDofGlobalApplication < this->m_dofBD[0/*iBD*/].numGlobalDofInterface(); ++idDofGlobalApplication, ++banner) {
116
1/2
✓ Branch 1 taken 324 times.
✗ Branch 2 not taken.
 324 this->clearMatrixRHS( FlagMatrixRHS::only_rhs );
117 // retrieving the correct volume indeces
118 324 idDofGlobalVolumePetsc=this->m_dofBD[0/*iBD*/].globBD2PetscVol(idDofGlobalApplication);
119 324 felInt idDofGlobalVolumeApplication = idDofGlobalVolumePetsc;
120 324 felInt idDofGlobalPetsc = idDofGlobalApplication;
121
1/2
✓ Branch 1 taken 324 times.
✗ Branch 2 not taken.
 324 AOPetscToApplication(this->m_ao,1,&idDofGlobalVolumeApplication);
122
1/2
✓ Branch 3 taken 324 times.
✗ Branch 4 not taken.
 324 AOApplicationToPetsc(this->m_dofBD[0/*iBD*/].ao(), 1, &idDofGlobalPetsc);
123
124 // creating the correct data steklov
125
3/6
✓ Branch 2 taken 324 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 324 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 324 times.
✗ Branch 9 not taken.
 324 this->m_vecs.Get("dataSteklov").set(0.0);
126
3/4
✓ Branch 2 taken 324 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✓ Branch 5 taken 243 times.
 324 if ( this->m_dofPart[idDofGlobalVolumeApplication] == MpiInfo::rankProc() ) {
127
3/6
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 81 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 9 not taken.
 81 this->m_vecs.Get("dataSteklov").setValue(idDofGlobalVolumePetsc, /* value */ 1., INSERT_VALUES);
128 }
129
3/6
✓ Branch 2 taken 324 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 324 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 324 times.
✗ Branch 9 not taken.
 324 this->m_vecs.Get("dataSteklov").assembly();
130
5/10
✓ Branch 2 taken 324 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 324 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 324 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 324 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 324 times.
✗ Branch 16 not taken.
 324 this->gatherVector(this->m_vecs.Get("dataSteklov"),this->m_seqVecs.Get("dataSteklov"));
131
132 //Problem solution
133
1/2
✓ Branch 1 taken 324 times.
✗ Branch 2 not taken.
 324 this->assembleVolumeSystem( FlagMatrixRHS::only_rhs );
134
135 324 m_chronoRS->start();
136
3/6
✓ Branch 1 taken 324 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 324 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 324 times.
✗ Branch 8 not taken.
 324 this->solve( MpiInfo::rankProc(), MpiInfo::numProc() );
137
1/2
✓ Branch 2 taken 324 times.
✗ Branch 3 not taken.
 324 m_chronoRS->stop();
138
139
1/2
✓ Branch 2 taken 324 times.
✗ Branch 3 not taken.
 324 if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
140
1/2
✓ Branch 7 taken 324 times.
✗ Branch 8 not taken.
 324 this->solution().getValues(this->m_dofBD[0/*iBD*/].numLocalDofInterface(), this->m_dofBD[0/*iBD*/].loc2PetscVolPtr(), m_tmpValues.data());
141
1/2
✓ Branch 6 taken 324 times.
✗ Branch 7 not taken.
 324 m_fullSteklov.setValues(/*nb rows*/this->m_dofBD[0/*iBD*/].numLocalDofInterface(),/*id rows*/this->m_dofBD[0/*iBD*/].loc2PetscBDPtr(), /*nb cols*/ 1,/*id cols*/ &idDofGlobalPetsc , m_tmpValues.data(), INSERT_VALUES);
142 }
143 }
144
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
145
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_fullSteklov.assembly(MAT_FINAL_ASSEMBLY);
146 }
147
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 banner.finalizeComputation();
148
149 // restoring the normal state
150
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 FelisceParam::linearSolverVerboseON();
151
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 this->useSteklovDataEnd();
152 }
153 4 break;
154 ✗ default:
155 ✗ FEL_ERROR("unrecognized option for Steklov operator");
156 }
157
158
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
159
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 4 switch ( FelisceParam::instance().optionForSteklov ) {
160 ✗ case LOAD_IT:
161 {
162 ✗ std::stringstream filename;
163 ✗ if ( FelisceParam::instance().steklovMatrixName.size() == 0 )
164 ✗ filename << "steklov_" << MpiInfo::numProc();
165 else
166 ✗ filename << FelisceParam::instance().steklovMatrixName << "_" << MpiInfo::numProc();
167 ✗ m_fullSteklov.loadFromBinaryFormat(this->m_dofBD[0/*iBD*/].comm(),filename.str(),FelisceParam::instance().steklovDataDir);
168 }
169 ✗ break;
170 4 case COMPUTE_IT:
171 4 break;
172 ✗ case COMPUTE_AND_SAVE_IT:
173 {
174 ✗ std::stringstream filename;
175 ✗ filename << "steklov_" << MpiInfo::numProc();
176 ✗ m_fullSteklov.saveInBinaryFormat(this->m_dofBD[0/*iBD*/].comm(),filename.str(),FelisceParam::instance().steklovDataDir);
177 }
178 ✗ break;
179 }
180 }
181
2/6
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 4 times.
 4 if ( FelisceParam::instance().exportSteklovMatrix && this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
182 ✗ m_fullSteklov.saveInBinaryFormat(this->m_dofBD[0/*iBD*/].comm(),"Steklov",FelisceParam::instance().resultDir);
183 }
184 4 }
185 /*! \brief Computation of the Reduced Order (Ro) Steklov operator.
186 *
187 * The function provides an approximation of the Steklov operator. The matrix rapresentation of the operator
188 * is never computated. The action of the operator on a new input std::vector is approximated via a low rank expansion.
189 */
190 template <class volumeProblem>
191 void
192 ✗ LinearProblemReducedSteklov<volumeProblem>::assembleLowRankSteklov(felInt imesh){
193 ✗ this->initChronoRS();
194 //! First the mass and the laplacian matrix on the boundary are computed.
195 ✗ PetscMatrix mass = assembleMassBoundary();
196 ✗ PetscMatrix laplacian = assembleLaplacianBoundary();
197 //! Then the object handling the steklov reduction techniques is created.
198 ✗ m_reducedSteklov = new ReducedSteklov<volumeProblem>(mass, laplacian, this->m_dofBD[0/*iBD*/].comm(), this, m_imgType,m_chronoRS);
199
200 ✗ this->computeTheConstantResponse();
201 ✗ m_reducedSteklov->addConstantResponse(this->getSteklovImg());
202
203 ✗ if ( m_reducedSteklov->loadedFromFile() ) {
204 ✗ std::stringstream msg;
205 ✗ msg<<"Off-line basis correctly loaded."<<std::endl;
206 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
207 ✗ } else {
208 ✗ if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
209 //! The laplacian eigen-problem on the surface is solved.
210 ✗ m_reducedSteklov->solveLaplacianEP();
211 }
212 //! The image through the steklov operator of the laplacian eigen-std::vector is computed.
213 ✗ m_reducedSteklov->applySteklovOperatorOnLEV(imesh);
214 ✗ if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
215 //! The steklov eigenValues and eigenVectors are approximated.
216 ✗ m_reducedSteklov->steklovEVecEstimate();
217 }
218 }
219 }
220 template <class volumeProblem>
221 void
222 ✗ LinearProblemReducedSteklov<volumeProblem>::setSteklovData( PetscVector& bdInput ) {
223 ✗ this->m_vecs.Get("dataSteklov").set(0.0);
224 ✗ this->dofBD(/*iBD*/0).extendOnVolume( this->m_vecs.Get("dataSteklov"), bdInput );
225 ✗ this->gatherVector( this->m_vecs.Get("dataSteklov"), this->m_seqVecs.Get("dataSteklov"));
226 }
227
228 template <class volumeProblem>
229 ✗ PetscVector LinearProblemReducedSteklov<volumeProblem>::getSteklovImg() {
230
231 //! Initialization
232 ✗ PetscVector bdOutput;
233 ✗ if ( m_tmpValues.size() != this->m_dofBD[0/*iBD*/].numLocalDofInterface() ) {
234 ✗ m_tmpValues.resize( this->m_dofBD[0/*iBD*/].numLocalDofInterface() );
235 }
236
237 ✗ if ( m_imgType == dirichlet ) {
238 //! Just read the data
239 ✗ this->solution().getValues( this->m_dofBD[0/*iBD*/].numLocalDofInterface(), this->m_dofBD[0/*iBD*/].loc2PetscVolPtr(), m_tmpValues.data());
240 ✗ } else if (m_imgType == neumann) {
241 //! Compute the residual
242 ✗ this->computeResidual();
243 //! And read the data inside
244 ✗ this->m_residual.scale(-1); //We change the sign to be consistent with the reduce steklov class
245 ✗ this->m_residual.getValues( this->m_dofBD[0/*iBD*/].numLocalDofInterface(), this->m_dofBD[0/*iBD*/].loc2PetscVolPtr(), m_tmpValues.data());
246 ✗ this->m_residual.scale(-1);
247 }
248 ✗ if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
249 //allocation
250 ✗ bdOutput = this->dofBD(/*iBD*/0).allocateBoundaryVector(DofBoundary::parallel);
251 //std::set to zero
252 ✗ bdOutput.set(0.0);
253 //std::set of the real values
254 ✗ bdOutput.setValues( this->m_dofBD[0/*iBD*/].numLocalDofInterface(), this->m_dofBD[0/*iBD*/].loc2PetscBDPtr(), m_tmpValues.data(), INSERT_VALUES);
255 //assembly
256 ✗ bdOutput.assembly();
257 }
258 ✗ return bdOutput;
259 }
260
261 // ======================================================================
262 // now the privates methods
263 // ======================================================================
264 // reader of vectors coming from different linearProblems
265 // the std::vector are stored in sequential std::vector "whereToSave"
266 // For scalar
267 template <class volumeProblem>
268 void
269 312 LinearProblemReducedSteklov<volumeProblem>::readerInterface( const std::map<int,std::map<int,int> > & externalMap,
270 PetscVector& seqExternalVec,
271 const std::map<int,int> &mapBetweenLabels,
272 std::string whereToSave,
273 felInt iUnkWTW, // std::set to zero by default
274 felInt iCompWTW // std::set to zero by default
275 ) {
276
1/2
✓ Branch 2 taken 312 times.
✗ Branch 3 not taken.
 312 this->readerInterfaceComponentWise(externalMap,seqExternalVec,mapBetweenLabels,whereToSave,iUnkWTW,iCompWTW);
277 312 this->m_vecs[whereToSave].assembly();
278 312 this->gatherVector(this->m_vecs[whereToSave],this->m_seqVecs[whereToSave]);
279
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 312 times.
 312 if ( FelisceParam::verbose() > 0 ) {
280 ✗ PetscPrintf(MpiInfo::petscComm(), "done\n");
281 }
282 312 }
283 // For std::vector
284 template <class volumeProblem>
285 void
286 ✗ LinearProblemReducedSteklov<volumeProblem>::readerInterface( const std::map<int, std::map<int,std::map<int,int> > > & externalMap,
287 PetscVector& seqExternalVec,
288 const std::map<int,int> &mapBetweenLabels,
289 std::string whereToSave,
290 felInt iUnkWTW,
291 felInt iUCR
292 ) {
293 ✗ for ( std::size_t iComp(0); iComp < (std::size_t) this->dimension(); ++iComp,++iUCR) {
294 ✗ this->readerInterfaceComponentWise( externalMap.at(iUCR),seqExternalVec, mapBetweenLabels, whereToSave, iUnkWTW, iComp);
295 }
296 ✗ this->m_vecs[whereToSave].assembly();
297 ✗ this->gatherVector(this->m_vecs[whereToSave],this->m_seqVecs[whereToSave]);
298 ✗ if ( FelisceParam::verbose() > 0 ) {
299 ✗ PetscPrintf(MpiInfo::petscComm(), "done\n");
300 }
301 }
302 // Reader of vectors coming from different linearProblems
303 // the std::vector are stored in sequential std::vector "whereToSave"
304 // TODO: This function fill a volume std::vector, can this be avoided?
305 template <class volumeProblem>
306 void
307 312 LinearProblemReducedSteklov<volumeProblem>::readerInterfaceComponentWise( const std::map<int,std::map<int,int> > & externalMap,
308 PetscVector& seqExternalVec,
309 const std::map<int,int> &mapBetweenLabels,
310 std::string whereToSave,
311 felInt iUnkWTW,
312 felInt iCompWTW) {
313 // some couts to help debugging
314
2/4
✓ Branch 1 taken 312 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 312 times.
 312 if ( FelisceParam::verbose() > 2 )
315 ✗ PetscPrintf(MpiInfo::petscComm(),"%s", std::string("[linearProblem--Reduced Steklov] Reading "+ whereToSave +"...").c_str() );
316
317 felInt iDofGlobThisLPB; // petsc index of the current dof in this linear problem
318 felInt iDofGlobThisLPBApp; // application index of the current dof in this linear problem
319 felInt iDofGlobOtherLPB; // petsc index of the current dof in the other linear problem
320
321
1/2
✓ Branch 2 taken 312 times.
✗ Branch 3 not taken.
 312 int iUnknComp = this->dof().getNumGlobComp(iUnkWTW, iCompWTW);
322 double tmpValue;
323
324 // for each label of the interface in the mesh of this problem
325
2/2
✓ Branch 4 taken 312 times.
✓ Branch 5 taken 312 times.
 624 for(auto currentLabel = m_interfaceLabels.begin(); currentLabel != m_interfaceLabels.end(); currentLabel++) {
326 // we get the corresponding label in the other mesh
327
1/2
✓ Branch 2 taken 312 times.
✗ Branch 3 not taken.
 312 felInt labelOtherLPBMesh=mapBetweenLabels.find(*currentLabel)->second;
328 // some couts to help debugging
329
2/4
✓ Branch 1 taken 312 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 312 times.
 312 if ( FelisceParam::verbose() > 2 ) {
330 ✗ std::cout << std::endl;
331 ✗ std::cout << "Interface label from reduced steklov mesh: "<<*currentLabel<<" , label from the other mesh: "<<labelOtherLPBMesh<<std::endl;
332 }
333 // for each dof on this portion of the interface
334
1/2
✓ Branch 1 taken 312 times.
✗ Branch 2 not taken.
 312 for(auto currentPetscOtherLPBDOF = externalMap.find(labelOtherLPBMesh)->second.begin(); //begin of the std::unordered_map
335
3/4
✓ Branch 1 taken 25584 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 25272 times.
✓ Branch 7 taken 312 times.
 25584 currentPetscOtherLPBDOF != externalMap.find(labelOtherLPBMesh)->second.end(); //end of the std::unordered_map
336 25272 ++currentPetscOtherLPBDOF) {
337 25272 iDofGlobOtherLPB = currentPetscOtherLPBDOF->second; //global dof
338
1/2
✓ Branch 4 taken 25272 times.
✗ Branch 5 not taken.
 25272 iDofGlobThisLPB = this->m_dofBD[0/*iBD*/].getPetscDofs(/*iUnknown,iComp where to write*/iUnknComp, *currentLabel, currentPetscOtherLPBDOF->first) ;
339
340 25272 iDofGlobThisLPBApp = iDofGlobThisLPB;
341
1/2
✓ Branch 1 taken 25272 times.
✗ Branch 2 not taken.
 25272 AOPetscToApplication(this->m_ao,1,&iDofGlobThisLPBApp);
342 // we check if the dof belongs to this proc
343
3/4
✓ Branch 2 taken 25272 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6318 times.
✓ Branch 5 taken 18954 times.
 25272 if ( this->m_dofPart[iDofGlobThisLPBApp] == MpiInfo::rankProc() ) {
344 //we get the value
345
1/2
✓ Branch 1 taken 6318 times.
✗ Branch 2 not taken.
 6318 seqExternalVec.getValues(1,&iDofGlobOtherLPB, &tmpValue);
346 // and we std::set it
347
3/6
✓ Branch 1 taken 6318 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6318 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6318 times.
✗ Branch 8 not taken.
 6318 this->m_vecs.Get(whereToSave).setValue(iDofGlobThisLPB, tmpValue, INSERT_VALUES);
348 }
349 }
350 }
351 312 }
352
353
354
355
356 template <class volumeProblem>
357 void
358 ✗ LinearProblemReducedSteklov<volumeProblem>::setValueMatrixBD(felInt ielSupportDof) {
359 // Elementary matrices are placed into the bigger one defined on the boundary.
360 ✗ felInt numDofTotal = this->m_elementMatBD[0]->mat().size1();
361
362 ✗ felInt node1 = 0;
363 ✗ felInt node2 = 0;
364 ✗ felInt cptGposLine = 0;
365 ✗ felInt cptGposLine2 = 0;
366 ✗ felInt sizeSupport = 0;
367 ✗ felInt sizeSupport2 = 0;
368 ✗ felInt cptBlock = 0;
369 ✗ int idVar = -1;
370 ✗ int idVar2 = -1;
371 //! for all the unknowns of the problems
372 ✗ for ( std::size_t iUnknown = 0; iUnknown < this->m_listUnknown.size(); iUnknown++) {
373 ✗ idVar = this->m_listUnknown.idVariable(iUnknown);
374 ✗ sizeSupport = this->m_listCurvilinearFiniteElement[idVar]->numDof();
375 //! for all its components
376 ✗ for (std::size_t iComp = 0; iComp < this->m_listVariable[idVar].numComponent(); iComp++) {
377 ✗ cptBlock++;
378 //! for all the dof on the element
379 ✗ for ( int iSupport = 0; iSupport < sizeSupport; iSupport++) {
380
381 //! get the volume application index
382 ✗ this->dof().loc2glob(ielSupportDof, iSupport, idVar, iComp, node1);
383 //! save it
384 ✗ m_globPosRow[cptGposLine]= node1;
385 ✗ cptGposLine2 = 0;
386 //! then we do the same with the columns:
387 //! for each unknown
388 ✗ for ( std::size_t iUnknown2 = 0; iUnknown2 <this->m_listUnknown.size(); iUnknown2++) {
389 ✗ idVar2 = this->m_listUnknown.idVariable(iUnknown2);
390 ✗ sizeSupport2 = this->m_listCurvilinearFiniteElement[idVar2]->numDof();
391 //! for all its components
392 ✗ for (std::size_t jComp = 0; jComp <this->m_listVariable[idVar2].numComponent(); jComp++) {
393 //! for all the dofs on the element
394 ✗ for ( int jSupport = 0; jSupport < sizeSupport2; jSupport++) {
395 ✗ int iConnect = this->dof().getNumGlobComp( iUnknown, iComp);
396 ✗ int jConnect = this->dof().getNumGlobComp( iUnknown2, jComp);
397 ✗ if ( this->m_listUnknown.mask()(iConnect, jConnect) > 0) {
398 //! the global index
399 ✗ this->dof().loc2glob(ielSupportDof,jSupport,idVar2,jComp,node2);
400 //! save it
401 ✗ m_globPosColumn[cptGposLine2] = node2;
402 //! store in this small matrix the values just computed
403 ✗ m_matrixValues[cptGposLine2 + cptGposLine*numDofTotal] = this->m_elementMatBD[0]->mat()(cptGposLine, cptGposLine2);
404 ✗ cptGposLine2++;
405 }
406 }
407 }
408 }
409 ✗ cptGposLine++;
410 }
411 }
412 }
413 //! Now we have to compute where this numbers have to be put into the boundary matrix
414 //! volume application -> petsc volume
415 ✗ AOApplicationToPetsc(this->m_ao,cptGposLine,m_globPosRow.data());
416 ✗ AOApplicationToPetsc(this->m_ao,cptGposLine2,m_globPosColumn.data());
417
418 ✗ std::vector<felInt> row;
419 ✗ std::vector<felInt> col;
420 ✗ std::vector<felInt> tmpCol;
421 ✗ std::vector<double> val;
422 //! petsc volume -> boundary application
423 ✗ for ( int k(0); k < cptGposLine; ++k) {
424 ✗ tmpCol.clear();
425 try {
426 // petscVol2ApplicationBD can throw an out_of_range while m_globPosRow[k] can not
427 // The idea is that if this fails it means that this particular dof is not in the mapping
428 // and that it can be discarded (e.g. pressure nodes when dealing with Stokes problem)
429 ✗ row.push_back(this->m_dofBD[0/*iBD*/].petscVol2ApplicationBD(m_globPosRow[k]));
430 ✗ for ( int k2(0); k2 < cptGposLine2; ++k2) {
431 try {
432 ✗ tmpCol.push_back(this->m_dofBD[0/*iBD*/].petscVol2ApplicationBD(m_globPosColumn[k2]));
433 ✗ val.push_back(m_matrixValues[k2+k*numDofTotal]);
434 ✗ } catch ( const std::out_of_range& oor ){}
435 }
436 ✗ } catch ( const std::out_of_range& oor ){}
437 ✗ if (tmpCol.size()>0) {
438 ✗ col=tmpCol;
439 }
440 }
441 //! boundary application -> petsc application
442 ✗ AOApplicationToPetsc(this->m_dofBD[0/*iBD*/].ao(),row.size(),row.data());
443 ✗ AOApplicationToPetsc(this->m_dofBD[0/*iBD*/].ao(),col.size(),col.data());
444 // std::set value
445 ✗ m_aux.setValues(row.size(),row.data(),col.size(),col.data(),val.data(),ADD_VALUES);
446 }
447
448
449 //! \brief It assembles the laplacian matrix on the boundary.
450 template <class volumeProblem>
451 PetscMatrix
452 ✗ LinearProblemReducedSteklov<volumeProblem>::assembleLaplacianBoundary() {
453 ✗ PetscMatrix laplacian;
454 ✗ if ( this->dimension() == 3 ) {
455 ✗ this->m_dofBD[0/*iBD*/].allocateMatrixOnBoundary( laplacian );
456 ✗ m_aux=laplacian;
457 ✗ if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
458 ✗ this->assemblyLoopBoundaryGeneral( &LinearProblemReducedSteklov<volumeProblem>::laplacianMatrixComputer,
459 ✗ m_interfaceLabels,
460 &LinearProblemReducedSteklov<volumeProblem>::initPerETLAP,
461 &LinearProblemReducedSteklov<volumeProblem>::updateFE);// we do not need the derivatives for surface laplacian: they are computed in m_curv
462 ✗ m_globPosColumn.clear();m_globPosRow.clear();m_matrixValues.clear();
463 ✗ laplacian.assembly(MAT_FINAL_ASSEMBLY);
464 }
465 ✗ if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
466 ✗ if ( FelisceParam::instance().exportLaplacianMatrix ) {
467 ✗ laplacian.saveInBinaryFormat( this->m_dofBD[0/*iBD*/].comm(), "Laplacian" , FelisceParam::instance().resultDir);
468 }
469 }
470 // Better to put it since m_aux will stay alive as long as the linear problem,
471 // while we would like to deallocate this matrix when the laplacian matrix gets out of scope
472 ✗ m_aux.destroy();
473 } else {
474 ✗ FEL_ERROR("We can not yet assemble the laplacian operator on the boundary of a 2D-domain.");
475 }
476 ✗ return laplacian;
477 }
478 /* ! \brief this function assemble the mass matrix at the boundary.
479 * It saves such matrix in binary format and then it destroies the matrix.
480 */
481 template <class volumeProblem>
482 PetscMatrix
483 ✗ LinearProblemReducedSteklov<volumeProblem>::assembleMassBoundary( ) {
484 ✗ PetscMatrix mass;
485 ✗ this->m_dofBD[0/*iBD*/].allocateMatrixOnBoundary( mass );
486 ✗ m_aux=mass;
487 ✗ if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
488 ✗ this->assemblyLoopBoundaryGeneral( &LinearProblemReducedSteklov<volumeProblem>::massMatrixComputer,
489 ✗ m_interfaceLabels,
490 &LinearProblemReducedSteklov<volumeProblem>::initPerETMASS,
491 &LinearProblemReducedSteklov<volumeProblem>::updateFE);
492 ✗ m_globPosColumn.clear();m_globPosRow.clear();m_matrixValues.clear();
493
494 ✗ mass.assembly(MAT_FINAL_ASSEMBLY);
495 ✗ if ( FelisceParam::instance().exportMassMatrix ) {
496 ✗ mass.saveInBinaryFormat( this->m_dofBD[0/*iBD*/].comm(), "Mass", FelisceParam::instance().resultDir);
497 }
498 }
499 // Better to put it since m_aux will stay alive as long as the linear problem,
500 // while we would like to deallocate this matrix when the mass matrix gets out of scope
501 ✗ m_aux.destroy();
502 ✗ return mass;
503 }
504 template <class volumeProblem>
505 void
506 ✗ LinearProblemReducedSteklov<volumeProblem>::exportOutputSteklov(FelisceTransient::Pointer fstransient,int iFixedPointIteration) {
507 ✗ PetscVector restrictedSolution = this->dofBD(/*iBD*/0).allocateBoundaryVector(DofBoundary::parallel);
508 ✗ this->dofBD(/*iBD*/0).restrictOnBoundary(this->solution(),restrictedSolution);
509 ✗ std::stringstream iFixedPointIterationString;
510 ✗ if ( iFixedPointIteration < 10 ) iFixedPointIterationString<<"00"<<iFixedPointIteration;
511 ✗ else if (iFixedPointIteration < 100) iFixedPointIterationString<<"0" <<iFixedPointIteration;
512 ✗ else iFixedPointIterationString <<iFixedPointIteration;
513 ✗ if( FelisceParam::instance().exportOutputOfSteklov ) {
514 // binary export
515 ✗ if ( FelisceParam::instance().exportOnlyLastOutput ) {
516 ✗ iFixedPointIterationString.str("last");
517 }
518 ✗ std::stringstream filename2;
519 ✗ filename2<<"outputSteklov_"<<fstransient->itStr()<<"_"<<iFixedPointIterationString.str().c_str();
520 ✗ if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
521 ✗ restrictedSolution.saveInBinaryFormat(this->m_dofBD[0/*iBD*/].comm(), filename2.str(), FelisceParam::instance().resultDir);
522 }
523 }
524 }
525
526 // the data is given as a volume std::vector: first you restrict this volume std::vector on
527 // the surface you apply the Steklov operator and you save the result in a volume
528 // std::vector
529 template <class volumeProblem>
530 void
531 304 LinearProblemReducedSteklov<volumeProblem>::applySteklov( PetscVector& volumeInput, PetscVector& volumeOutput, FelisceTransient::Pointer fstransient, int iFixedPointIteration) {
532 // reset this to zero
533
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 304 volumeOutput.set(0.0);
534
1/2
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
 304 PetscVector Input = this->dofBD(/*iBD*/0).allocateBoundaryVector(DofBoundary::parallel);
535
1/2
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
 304 PetscVector Output = this->dofBD(/*iBD*/0).allocateBoundaryVector(DofBoundary::parallel);
536
1/2
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
 304 this->dofBD(/*iBD*/0).restrictOnBoundary( volumeInput, Input);
537
538
1/2
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
 304 if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) {
539
2/4
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 152 times.
 304 if ( FelisceParam::instance().useRoSteklov ) {
540 ✗ m_reducedSteklov->applyLowRankSteklov(Input,Output);
541 } else {
542
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 304 mult(m_fullSteklov, Input, Output);
543 }
544
1/2
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
 304 this->exportSteklovInputOutput(Input,Output,fstransient,iFixedPointIteration);
545 }
546
1/2
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
 304 this->dofBD(/*iBD*/0).extendOnVolume(volumeOutput,Output);
547 }
548 template <class volumeProblem>
549 304 void LinearProblemReducedSteklov<volumeProblem>::exportSteklovInputOutput( PetscVector& Input, PetscVector& Output, FelisceTransient::Pointer fstransient, int iFixedPointIteration) {
550
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 152 times.
 304 if ( FelisceParam::instance().exportSteklovData ) {
551 ✗ std::stringstream iFixedPointIterationString;
552 ✗ if ( iFixedPointIteration < 10 ) iFixedPointIterationString<<"00"<<iFixedPointIteration;
553 ✗ else if (iFixedPointIteration < 100) iFixedPointIterationString<<"0" <<iFixedPointIteration;
554 ✗ else iFixedPointIterationString <<iFixedPointIteration;
555 ✗ if( FelisceParam::instance().exportInputOfSteklov ) {
556 ✗ if ( FelisceParam::instance().exportOnlyLastInput ) { iFixedPointIterationString.str("last"); }
557 ✗ std::stringstream filename;
558 ✗ filename<<"inputSteklov_"<<fstransient->itStr()<<"_"<<iFixedPointIterationString.str().c_str();
559 ✗ Input.saveInBinaryFormat(this->m_dofBD[0/*iBD*/].comm(), filename.str(), FelisceParam::instance().resultDir);
560 }
561 ✗ if( FelisceParam::instance().exportOutputOfSteklov ) {
562 ✗ if ( FelisceParam::instance().exportOnlyLastOutput ) { iFixedPointIterationString.str("last"); }
563 ✗ std::stringstream filename2;
564 ✗ filename2<<"outputSteklov_"<<fstransient->itStr()<<"_"<<iFixedPointIterationString.str().c_str();
565 ✗ Output.saveInBinaryFormat(this->m_dofBD[0/*iBD*/].comm(), filename2.str(), FelisceParam::instance().resultDir);
566 }
567 }
568 }
569 template <class volumeProblem>
570 ✗ void LinearProblemReducedSteklov<volumeProblem>::userSetNullSpace(PetscVector& v, int k) {
571 //TODO move this in a user file or do something better
572 ✗ if ( FelisceParam::instance().notConnectedCylinders ) {
573 ✗ v.set(0.0);
574 ✗ std::size_t labelCount = 0;
575 ✗ double aus(1.0);
576 ✗ for(auto itLabel = this->m_dofBD[0/*iBD*/].listOfBoundaryPetscDofs(/*iUnkown*/0,/*iComp*/0).begin(); itLabel != this->m_dofBD[0/*iBD*/].listOfBoundaryPetscDofs(/*iUnkown*/0,/*iComp*/0).end(); ++itLabel, ++labelCount ) {
577 ✗ if ( labelCount == std::size_t(k) ) {
578 ✗ for(auto itDof = itLabel->second.begin(); itDof != itLabel->second.end(); ++itDof) {
579 ✗ felInt dof=this->m_dofBD[0/*iBD*/].petscVol2ApplicationBD(itDof->second);
580 ✗ AOApplicationToPetsc(this->m_dofBD[0/*iBD*/].ao(),1,&dof);
581 ✗ v.setValues(1,&dof,&aus, INSERT_VALUES);
582 }
583 }
584 }
585 ✗ v.assembly();
586 } else {
587 ✗ if ( k > 0 ) {
588 ✗ FEL_ERROR("There are too many zero eigenvalues: implement your own userSetNullSpace function");
589 } else {
590 ✗ v.set(1.0);
591 }
592 }
593 }
594 //! \brief it modifies the state of the class
595 //! In particular it raises the flag m_useSteklovData, and it allocates the dataSteklov vectors
596 template <class volumeProblem>
597 4 void LinearProblemReducedSteklov<volumeProblem>::useSteklovDataBegin() {
598 4 this->m_useSteklovData = true;
599
3/6
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
 4 this->m_vecs["dataSteklov"].duplicateFrom( this->solution() );
600
3/6
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
 4 this->m_vecs.Get("dataSteklov").set(0.0);
601
3/6
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
 4 this->m_seqVecs["dataSteklov"].duplicateFrom( this->sequentialSolution() );
602
3/6
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
 4 this->m_seqVecs.Get("dataSteklov").set(0.0);
603 4 m_tmpValues.resize( this->m_dofBD[0/*iBD*/].numLocalDofInterface() );
604 4 this->m_boundaryConditionList.temporarySetBCValuesToZero();
605 4 }
606
607 //! \brief it restores the previous state of the class
608 //! Call this function only after having used useSteklovBegin.
609 //! The flag: m_useSteklovData is reset to false and the dataSteklov vectors are destroyed
610 template <class volumeProblem>
611 void
612 4 LinearProblemReducedSteklov<volumeProblem>::useSteklovDataEnd() {
613 4 this->m_useSteklovData = false;
614
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 this->m_vecs.erase("dataSteklov");
615
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 this->m_seqVecs.erase("dataSteklov");
616 4 this->m_boundaryConditionList.restoreBCValues();
617 4 }
618 }
619