GCC Code Coverage Report


Directory: ./
File: PETScInterface/romPETSC.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 231 0.0%
Branches: 0 410 0.0%

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: E. Schenone
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "PETScInterface/romPETSC.hpp"
21
22 namespace felisce {
23
24 RomPETSC::RomPETSC():
25 m_fstransient(nullptr),
26 m_buildBasis(false),
27 m_buildMass(false),
28 m_buildInverse(false),
29 m_isFirstStep(false)
30 {
31 m_setPreconditionerOption = SAME_NONZERO_PATTERN;
32 }
33
34 RomPETSC::~RomPETSC() {
35 if(m_buildBasis) {
36 m_romBasis.destroy();
37 m_romMatrix.destroy();
38 }
39 if(m_buildMass) {
40 m_mass.destroy();
41 }
42 if(m_buildInverse) {
43 m_inverseMatrix.destroy();
44 }
45 if (m_isFirstStep) {
46 m_romRHS.destroy();
47 m_romSolution.destroy();
48 }
49 }
50
51 // initialize the pod basis parameters (passing AO& ao too - useful to find petsc dof indeces mapping)
52 void RomPETSC::initializeRom(FelisceTransient::Pointer fstransient,MPI_Comm& comm, AO& ao) {
53 m_fstransient = fstransient;
54 m_verbose = FelisceParam::verbose();
55 m_MPIcomm = comm;
56 m_ao = ao;
57 m_filenameBasis = FelisceParam::instance().podBasisFile;
58 m_dimRomBasis = FelisceParam::instance().dimRomBasis;
59 m_dimFullRomBasis = m_dimRomBasis;
60 m_useInverse = true;
61 }
62
63 // initialize the pod basis parameters
64 void RomPETSC::initializeRom(FelisceTransient::Pointer fstransient,MPI_Comm& comm) {
65 m_fstransient = fstransient;
66 m_verbose = FelisceParam::verbose();
67 m_MPIcomm = comm;
68 m_filenameBasis = FelisceParam::instance().podBasisFile;
69 m_dimRomBasis = FelisceParam::instance().dimRomBasis;
70 m_dimFullRomBasis = m_dimRomBasis;
71 m_useInverse = true;
72 }
73
74 // Intialize basis matrix from matrix
75 void RomPETSC::initializeRomBasis(const PetscMatrix& basis) {
76 felInt lengthBasis;
77 basis.getSize(&lengthBasis,&m_dimRomBasis);
78
79 m_romBasis.duplicateFrom(basis,MAT_COPY_VALUES);
80 m_romBasis.assembly(MAT_FINAL_ASSEMBLY);
81
82 m_buildBasis = true;
83
84 }
85
86 // Intialize basis matrix from matrix
87 void RomPETSC::initializeRomBasis(const std::vector<PetscVector*>& basis) {
88 m_dimRomBasis = basis.size();
89 felInt lengthBasis;
90 (*basis[0]).getSize(&lengthBasis);
91 m_romBasis.createDense(m_MPIcomm,PETSC_DECIDE,PETSC_DECIDE,lengthBasis,m_dimRomBasis,FELISCE_PETSC_NULLPTR);
92 m_romBasis.setFromOptions();
93
94 double tmpValue[lengthBasis];
95 felInt ia[lengthBasis];
96 for (felInt i=0; i< lengthBasis; i++) {
97 ia[i] = i;
98 }
99 AOApplicationToPetsc(m_ao,lengthBasis,ia);
100 int ja;
101 for (int j=0; j<m_dimRomBasis; j++) {
102 (*basis[j]).getValues(lengthBasis,ia,tmpValue);
103 ja = j;
104 AOApplicationToPetsc(m_ao,1,&ja);
105 m_romBasis.setValues(lengthBasis,ia,1,&ja,tmpValue,INSERT_VALUES);
106 }
107 m_romBasis.assembly(MAT_FINAL_ASSEMBLY);
108
109 m_buildBasis = true;
110
111 }
112
113 // read POD basis from file
114 // format:
115 // # of elements
116 // # of dof
117 // element1 || ...|| elementK
118 // element1 || ...|| elementK
119 // ...
120 // Intialize basis matrix from file
121 void RomPETSC::initializeRomBasis() {
122 m_filenameBasis = FelisceParam::instance().inputDirectory + "/" + FelisceParam::instance().podBasisFile;
123 if (m_verbose) {
124 std::cout << " Rom::initBasis: Reading from file "<< m_filenameBasis << std::endl;
125 }
126 std::ifstream basisFile(m_filenameBasis.c_str(),std::ios_base::in);
127 if (!basisFile) {
128 FEL_ERROR(" Rom::initBasis(): Impossible to read file: "+m_filenameBasis+" (POD basis)");
129 }
130
131 felInt lengthBasis;
132 basisFile >> lengthBasis; // read n. of rows
133 basisFile >> m_dimFullRomBasis; // read n. of columns (total)
134 if (m_verbose) {
135 std::cout << "Read matrix of " << lengthBasis << " rows and " << m_dimFullRomBasis << " columns.\n";
136 }
137
138 if (m_dimFullRomBasis<m_dimRomBasis) {
139 std::cout << " requested dimBasis = " << m_dimRomBasis << ", available elements = "<< m_dimFullRomBasis << std::endl;
140 FEL_ERROR(" Rom::initBasis: not enough basis element available !");
141 }
142
143 m_romBasis.createDense(m_MPIcomm, PETSC_DECIDE, PETSC_DECIDE, lengthBasis, m_dimRomBasis, FELISCE_PETSC_NULLPTR);
144 m_romBasis.setFromOptions();
145
146 double tmpValue[m_dimRomBasis];
147 felInt ia[m_dimRomBasis];
148 felInt ja;
149 char line[1024];
150 double k = 1.0;
151 for (felInt j=0; j<lengthBasis; j++) { //rows
152 ja = j;
153 if ((j+1.0)>=(k*0.1*lengthBasis)) {
154 std::cout << "... " << 10*k << "% ";
155 k=k+1.0;
156 }
157 for (felInt i=0; i< m_dimRomBasis; i++) { // columns
158 ia[i] = i;
159 basisFile >> tmpValue[i];
160 }
161 basisFile.getline(line, 1024, '\n');
162 AOApplicationToPetsc(m_ao,1,&ja);
163 m_romBasis.setValues(1,&ja,m_dimRomBasis,ia,tmpValue,INSERT_VALUES);
164 }
165 m_romBasis.assembly(MAT_FINAL_ASSEMBLY);
166
167 m_buildBasis = true;
168 }
169
170 void RomPETSC::gatherVector(PetscVector& parVector, PetscVector& seqVector) {
171 PetscVector vout;
172 parVector.scatterToAll(vout,INSERT_VALUES,SCATTER_FORWARD);
173 seqVector.duplicateFrom(vout);
174 seqVector.copyFrom(vout);
175 vout.destroy();
176 }
177
178 // given a full std::vector f, compute \hat f = Phi^t f
179 void RomPETSC::fullToReducedRHS(PetscVector& fullVector) {
180 if (m_fstransient->iteration == 1) {
181 m_romRHS.create(m_MPIcomm);
182 m_romRHS.setSizes(PETSC_DECIDE, m_dimRomBasis);
183 m_romRHS.setFromOptions();
184 m_romRHS.zeroEntries();
185 }
186
187 multTranspose(m_romBasis,fullVector,m_romRHS);
188
189 m_isFirstStep = true;
190 }
191
192 // given a full matrix A, compute \hat A = Phi^t A \Phi
193 void RomPETSC::fullToReducedMatrix(PetscMatrix& fullMatrix) {
194 PetscMatrix tmpMat;
195 // m_tmpMat = fullMatrix * m_romBasis
196 matMult(fullMatrix, m_romBasis, MAT_INITIAL_MATRIX, PETSC_DEFAULT, tmpMat);
197 // m_romMatrix = m_romBasis^T * tmpMat
198 transposeMatMult(m_romBasis, tmpMat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, m_romMatrix);
199
200 tmpMat.destroy();
201
202 if (m_useInverse)
203 calculateInverse();
204
205 }
206
207 // solve the reduced linear system
208 void RomPETSC::solveALP() {
209 if (m_fstransient->iteration == 1) {
210 m_romSolution.duplicateFrom(m_romRHS);
211 }
212 m_romSolution.zeroEntries();
213
214 if (m_verbose) {
215 std::cout << " Rom:: solveROM() - solve reduced model "<< std::endl;
216 }
217
218 mult(m_inverseMatrix, m_romRHS, m_romSolution);
219
220 }
221
222
223 // solve the reduced linear system
224 void RomPETSC::solveROM(PetscVector& solution) {
225 if (m_fstransient->iteration == 1) {
226 m_romSolution.duplicateFrom(m_romRHS);
227 }
228 m_romSolution.zeroEntries();
229
230 if (m_verbose) {
231 std::cout << " Rom:: solveROM() - solve reduced model "<< std::endl;
232 }
233
234 if (!m_useInverse) {
235
236 PCFactorSetReuseFill(m_pc, PETSC_TRUE);
237
238 KSPSetOperators(m_ksp,m_romMatrix.toPetsc(),m_romMatrix.toPetsc());
239 if(strcmp(FelisceParam::instance().setPreconditionerOption[0].c_str(), "SAME_PRECONDITIONER") == 0)
240 KSPSetReusePreconditioner(m_ksp, PETSC_TRUE);
241 else
242 KSPSetReusePreconditioner(m_ksp, PETSC_FALSE);
243
244 KSPSetFromOptions(m_ksp);
245 KSPSolve(m_ksp,m_romRHS.toPetsc(),m_romSolution.toPetsc());
246
247 PetscInt its;
248 KSPGetIterationNumber(m_ksp,&its);
249 if (m_verbose) {
250 PetscPrintf(m_MPIcomm, "\n/================Information about solver==============/\n");
251 PetscPrintf(m_MPIcomm, "Solver use: <%s>.\n",FelisceParam::instance().solver[0].c_str());
252 PetscPrintf(m_MPIcomm, "Preconditioner use: <%s>.\n",FelisceParam::instance().preconditioner[0].c_str());
253 PetscPrintf(m_MPIcomm, "Preconditioner re-use option: <%s>.\n",FelisceParam::instance().setPreconditionerOption[0].c_str());
254 if (m_verbose > 1 ) {
255 PetscPrintf(m_MPIcomm, "Relative tolerance: %lf.\n",FelisceParam::instance().relativeTolerance[0]);
256 PetscPrintf(m_MPIcomm, "Absolute tolerance: %lf.\n",FelisceParam::instance().absoluteTolerance[0]);
257 PetscPrintf(m_MPIcomm, "Maximum of iterations: %d.\n",FelisceParam::instance().maxIteration[0]);
258 }
259 PetscPrintf(m_MPIcomm, "Number of iterations to inverse matrix: %d\n", its);
260 }
261 } else {
262 mult(m_inverseMatrix, m_romRHS, m_romSolution);
263 }
264
265 if (m_verbose) {
266 std::cout << " Rom:: solveROM() - get full solution "<< std::endl;
267 }
268 mult(m_romBasis,m_romSolution,solution);
269 }
270
271 void RomPETSC::setSolver(KSP ksp, PC pc, MatStructure setPreconditionerOption)
272 {
273 IGNORE_UNUSED_ARGUMENT(ksp);
274 IGNORE_UNUSED_ARGUMENT(pc);
275 IGNORE_UNUSED_ARGUMENT(setPreconditionerOption);
276
277 const auto& r_instance = FelisceParam::instance();
278 KSPCreate(m_MPIcomm,&m_ksp);
279 PCCreate(m_MPIcomm,&m_pc);
280 const auto& solver = r_instance.solver[0];
281 if ( solver == "preonly" || solver == "mumps" || solver == "superlu_dist" )
282 KSPSetType(m_ksp, KSPPREONLY);
283 else
284 KSPSetType(m_ksp, solver.c_str());
285 PCSetFromOptions(m_pc);
286 KSPGetPC(m_ksp,&m_pc);
287 PCSetType(m_pc,r_instance.preconditioner[0].c_str());
288 KSPSetTolerances(m_ksp,r_instance.relativeTolerance[0],r_instance.absoluteTolerance[0],PETSC_DEFAULT,r_instance.maxIteration[0]);
289
290 // TODO: MATSOLVERPASTIX once it works
291 if ( solver == "preonly" || solver == "mumps" ) {
292 PCFactorSetMatSolverType(m_pc, MATSOLVERMUMPS);
293 } else
294 #if SUPERLU_DIST_COMPILED
295 if ( solver == "superlu_dist" ) {
296 PCFactorSetMatSolverType(m_pc, MATSOLVERSUPERLU_DIST);
297 } else
298 #else
299 if ( solver == "superlu_dist" ) {
300 FEL_WARNING("SuperLUDist not available, we proceed using MUMPS")
301 PCFactorSetMatSolverType(m_pc, MATSOLVERMUMPS);
302 } else
303 #endif
304 {
305 if ( r_instance.initSolverWithPreviousSolution[0] ) {
306 KSPSetInitialGuessNonzero(m_ksp, PETSC_TRUE);
307 }
308 }
309
310 {
311 // TODO: to be removed
312 PetscOptionsSetValue(NULL, "-mat_mumps_icntl_20", "0");
313 }
314 }
315
316 void RomPETSC::calculateInverse() {
317 //Create the identityMatrix
318 PetscMatrix identityMatrix;
319 identityMatrix.createDense(m_MPIcomm, PETSC_DECIDE, PETSC_DECIDE, m_dimRomBasis, m_dimRomBasis, FELISCE_PETSC_NULLPTR);
320 identityMatrix.setFromOptions();
321 identityMatrix.zeroEntries();
322 m_inverseMatrix.duplicateFrom(identityMatrix,MAT_COPY_VALUES);
323 identityMatrix.shift(1.0);
324
325 // ILU factorization of_romMatrix
326 IS perm;
327 ISCreateStride(m_MPIcomm, m_dimRomBasis, 0, 1, &perm);
328 MatFactorInfo iluinfo;
329 iluinfo.fill = 1.0;
330 m_romMatrix.luFactor(perm, perm, &iluinfo);
331
332 // Solve m_romMatrix * m_inverseMatrix = identityMatrix
333 MatMatSolve(m_romMatrix.toPetsc(), identityMatrix.toPetsc(), m_inverseMatrix.toPetsc());
334
335 identityMatrix.destroy();
336
337 m_buildInverse = true;
338
339 }
340
341
342 void RomPETSC::readCase(int& precisionSol,std::vector<std::string>& variable, felInt& snapSize) {
343 // Solution files.
344 std::string caseFileName = FelisceParam::instance().resultDir + "/" + FelisceParam::instance().prefixName[0] + ".case";
345
346 std::ifstream casefile(caseFileName.c_str(), std::ios_base::in);
347 char line[1024];
348 std::string aux;
349 if ( !casefile ) {
350 std::cout << " ERROR: Can not open file "+caseFileName+"."<< std::endl;
351 exit(1);
352 } else {
353 if (FelisceParam::verbose() > 0) {
354 std::cout << "Reading " << caseFileName << std::endl;
355 }
356 }
357
358 casefile.getline(line, 1024, ':');
359 casefile.getline(line, 1024, ':');
360 casefile.getline(line, 1024, ' ');
361 casefile.getline(line, 1024, ' ');
362 casefile.getline(line, 1024, '\n');
363 std::string geoFileName = FelisceParam::instance().resultDir + "/" + line;
364 std::ifstream geoFile(geoFileName.c_str(),std::ios_base::in);
365 for (int i=0; i<5; i++) {
366 geoFile.getline(line,1024,'\n');
367 }
368 geoFile >> snapSize;
369 geoFile.close();
370
371 bool endCase = false;
372 while (!endCase) {
373 casefile.getline(line, 1024, ':');
374 casefile.getline(line, 1024, ' ');
375 casefile.getline(line, 1024, ' ');
376 casefile.getline(line, 1024, ' ');
377 casefile.getline(line, 1024, '.');
378 variable.emplace_back(line);
379 casefile.getline(line, 1024, '.');
380 aux = line;
381 precisionSol = aux.length();
382 casefile.getline(line, 1024, '\n');
383 casefile.getline(line, 1024, ' ');
384 aux = line;
385 if (aux == "TIME\ntime") {
386 endCase = true;
387 }
388 }
389 casefile.close();
390
391 if (m_verbose) {
392 std::cout << "Variables:\n";
393 for (std::size_t k=0; k < variable.size(); k++) {
394 std::cout << variable[k] << std::endl;
395 }
396 std::cout << "Numbers of nodes per variable : " << snapSize << std::endl;
397 }
398
399 }
400
401 void RomPETSC::readEnsight(double* scalarValue,std::string sclFile, felInt snapSize) {
402 std::string snapFileName = FelisceParam::instance().resultDir + sclFile + ".scl";
403 std::ifstream snap(snapFileName.c_str(), std::ios_base::in);
404
405 if ( !snap ) {
406 std::cout << " ERROR: Can not open file "+snapFileName+"."<< std::endl;
407 } else {
408 if (FelisceParam::verbose() > 1) {
409 std::cout << "Reading " << snapFileName << std::endl;
410 }
411 }
412 char newline[1024];
413 snap.getline(newline, 1024, '\n');
414 for (felInt j=0; j < snapSize; j++) {
415 snap >> scalarValue[j];
416 }
417 snap.close();
418 }
419
420
421 }
422