GCC Code Coverage Report

Directory: ./
File: Solver/coefProblem.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 94 0.0%
Branches: 0 104 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: J Fouchet & C Bui
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Core/felisceTools.hpp"
21 #include "Solver/coefProblem.hpp"
22 #include "PETScInterface/petscVector.hpp"
23
24 namespace felisce {
25 /*! Construtor */
26 ✗ CoefProblem::CoefProblem():
27 ✗ m_verbose(0),
28 ✗ m_buildSolver(false)
29 {}
30
31 ✗ CoefProblem::~CoefProblem() {
32 ✗ if(m_buildSolver) {
33 }
34 ✗ _RHS.destroy();
35 ✗ m_sol.destroy();
36 ✗ m_Matrix.destroy();
37 ✗ timeStepOverC.clear();
38 ✗ RplusTimeStepOverC.clear();
39 }
40
41 ✗ void CoefProblem::initializeCoefProblem(MPI_Comm& comm) {
42 ✗ m_petscComm = comm;
43 ✗ m_verbose = FelisceParam::verbose();
44
45 ✗ timeStepOverC.resize(FelisceParam::instance().lumpedModelBC_Rprox.size(), 0.);
46 ✗ RplusTimeStepOverC.resize(FelisceParam::instance().lumpedModelBC_Rprox.size());
47 ✗ RplusTimeStepOverC = FelisceParam::instance().lumpedModelBC_Rprox;
48
49 ✗ for(unsigned int i = 0; i < FelisceParam::instance().lumpedModelBC_Rprox.size(); i++) {
50 ✗ if ( FelisceParam::instance().lumpedModelBCType[i] == 1 ) {
51 // R lumpedModel bc, timeStepOverC = 0, RplusTimeStepOverC = R
52 // RCR lumpedModel bc: not implemented yet
53 ✗ if ( Tools::notEqual(FelisceParam::instance().lumpedModelBC_Rdist[i],0.) )
54 ✗ FEL_WARNING("Implicit lumpedModelBC are only implemented for R and RC elements.\n");
55 ✗ } else if ( FelisceParam::instance().lumpedModelBCType[i] == 2 ) {
56 // RC lumpedModel bc, timeStepOverC = dt/C, RplusTimeStepOverC = R + timeStepOverC
57 ✗ timeStepOverC[i] = FelisceParam::instance().timeStep / FelisceParam::instance().lumpedModelBC_C[i];
58 ✗ RplusTimeStepOverC[i] = FelisceParam::instance().lumpedModelBC_Rprox[i] + timeStepOverC[i];
59 }
60 }
61 }
62
63 /*!
64 \brief use the pattern define with dof to allocate memory for the matrix in CSR format.
65 */
66 ✗ void CoefProblem::allocateMatrixRHS(int size, int rank) {
67 IGNORE_UNUSED_RANK;
68
69 // Matrix
70 //========
71
72 ✗ m_Matrix.createDense(m_petscComm, PETSC_DECIDE, PETSC_DECIDE, size, size, FELISCE_PETSC_NULLPTR);
73 ✗ m_Matrix.zeroEntries();
74
75 // RHS
76 //=====
77
78 ✗ _RHS.createMPI(m_petscComm,PETSC_DECIDE,size);
79 ✗ m_sol.createMPI(m_petscComm,PETSC_DECIDE,size);
80 ✗ m_sol.zeroEntries();
81 }
82
83 ✗ void CoefProblem::assembleMatrix(const int idInflowOutflowBC, const std::vector<double> m_fluxLinCombBC, int rank) {
84 IGNORE_UNUSED_RANK;
85 // Matrix (I-A) is filled by columns (a "u_i" by column)
86 // We fill the "idInflowOutflowBC"th column (corresponding to the "idInflowOutflowBC"th Neumann Normal boundary condition)
87
88 ✗ PetscInt numLumpedModelBC = FelisceParam::instance().lumpedModelBC_Rprox.size() ;
89 ✗ PetscInt numNeumannNormalBCwithoutWindkessel = m_fluxLinCombBC.size() - numLumpedModelBC;
90 PetscScalar value;
91
92 // loop on the row of this column
93 // 1/ NeumannNormal boundary conditions without windkessel one (no resistance at the inlet)
94 ✗ for(PetscInt idRow = 0; idRow < numNeumannNormalBCwithoutWindkessel; idRow++) {
95 ✗ if (idRow == idInflowOutflowBC)
96 ✗ value = 1.;
97 else
98 ✗ value = 0.;
99 ✗ m_Matrix.setValues( 1, &idRow, 1, &idInflowOutflowBC, &value, INSERT_VALUES);
100 }
101
102 // 2/ windkessel bc
103 ✗ for(PetscInt idRow = numNeumannNormalBCwithoutWindkessel; idRow < (PetscInt) m_fluxLinCombBC.size(); idRow++) {
104 ✗ if (idRow == idInflowOutflowBC)
105 ✗ value = 1. - RplusTimeStepOverC[idRow-numNeumannNormalBCwithoutWindkessel] * m_fluxLinCombBC[idRow];
106 else
107 ✗ value = - RplusTimeStepOverC[idRow-numNeumannNormalBCwithoutWindkessel] * m_fluxLinCombBC[idRow];
108 ✗ m_Matrix.setValues( 1, &idRow, 1, &idInflowOutflowBC, &value, INSERT_VALUES);
109 }
110
111 ✗ m_Matrix.assembly(MAT_FINAL_ASSEMBLY);
112 }
113
114 ✗ void CoefProblem::assembleRHS(const std::vector<double> pressure, const std::vector<double> m_fluxLinCombBC, const std::vector<double> m_fluxLinCombBCtotal, int rank) {
115 IGNORE_UNUSED_RANK;
116 ✗ PetscInt numLumpedModelBC = FelisceParam::instance().lumpedModelBC_Rprox.size() ;
117 ✗ PetscInt numNeumannNormalBCwithoutWindkessel = m_fluxLinCombBC.size() - numLumpedModelBC;
118 PetscScalar value;
119
120 // loop on the row of this column
121 // 1/ NeumannNormal boundary conditions without windkessel one (no resistance at the inlet)
122 ✗ for(PetscInt idRow = 0; idRow < numNeumannNormalBCwithoutWindkessel; idRow++) {
123 // warning : b[i] = P_i and not -P_i
124 ✗ value = - pressure[idRow];
125 ✗ _RHS.setValues(1, &idRow, &value, INSERT_VALUES);
126 }
127
128 // 2/ windkessel bc
129 ✗ for(PetscInt idRow = numNeumannNormalBCwithoutWindkessel; idRow < (PetscInt) m_fluxLinCombBC.size(); idRow++) {
130 ✗ value = RplusTimeStepOverC[idRow-numNeumannNormalBCwithoutWindkessel] * m_fluxLinCombBC[idRow] + timeStepOverC[idRow-numNeumannNormalBCwithoutWindkessel] * m_fluxLinCombBCtotal[idRow];
131 ✗ _RHS.setValues(1, &idRow, &value, INSERT_VALUES);
132 }
133
134 ✗ _RHS.assembly();
135 }
136
137 /*!
138 \brief solve the linear solver with Petsc methods.
139 */
140 ✗ void CoefProblem::buildSolver()
141 {
142 ✗ m_ksp->init();
143 ✗ m_ksp->setKSPandPCType("preonly", "lu");
144 ✗ m_ksp->setTolerances(10.e-6, 10.e-10, 1000, 0);
145 ✗ m_buildSolver = true;
146 }
147
148 ✗ void CoefProblem::solve( int /*rank*/, int /*size*/ )
149 {
150 ✗ m_ksp->setKSPOperator(m_Matrix, FelisceParam::instance().setPreconditionerOption[0]);
151 ✗ m_ksp->solve(_RHS,m_sol,0);
152
153 ✗ printSolution(m_verbose);
154 }
155
156 // Write
157 //==========
158
159 //Write in file matrix.m, matrix in matlab format.
160 ✗ void CoefProblem::writeMatrixForMatlab(int verbose, std::ostream& outstr) const {
161 IGNORE_UNUSED_OUTSTR;
162 ✗ if (verbose > 20) {
163 PetscViewer viewer;
164 ✗ PetscViewerCreate(m_petscComm,&viewer);
165 ✗ PetscViewerASCIIOpen(m_petscComm,"__matrix.m",&viewer);
166 ✗ PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
167 ✗ m_Matrix.view(viewer);
168 }
169 }
170
171 //Write in file RHS.m, RHS in matlab format.
172 ✗ void CoefProblem::writeRHSForMatlab(int verbose, std::ostream& outstr) const {
173 IGNORE_UNUSED_OUTSTR;
174 ✗ if (verbose > 20) {
175 PetscViewer viewer;
176 ✗ PetscViewerCreate(m_petscComm,&viewer);
177 ✗ PetscViewerASCIIOpen(m_petscComm,"__RHS.m",&viewer);
178 ✗ PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
179 ✗ _RHS.view(viewer);
180 }
181 }
182
183 ✗ void CoefProblem::writeSolForMatlab(int verbose, std::ostream& outstr) const {
184 IGNORE_UNUSED_OUTSTR;
185 ✗ if (verbose > 20) {
186 PetscViewer viewer;
187 ✗ PetscViewerCreate(m_petscComm,&viewer);
188 ✗ PetscViewerASCIIOpen(m_petscComm,"m_m_sol.m",&viewer);
189 ✗ PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
190 ✗ m_sol.view(viewer);
191 }
192 }
193
194 //Write in file filename, matrix in matlab format.
195 ✗ void CoefProblem::writeMatrixForMatlab(std::string const& fileName, PetscMatrix& matrix) const {
196 PetscViewer viewer;
197 ✗ PetscViewerCreate(m_petscComm,&viewer);
198 ✗ PetscViewerASCIIOpen(m_petscComm,fileName.c_str(),&viewer);
199 ✗ PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
200 ✗ matrix.view(viewer);
201 }
202
203 //Write in file filenam, RHS in matlab format.
204 ✗ void CoefProblem::writeVectorForMatlab(std::string const& fileName, PetscVector& vector) const {
205 PetscViewer viewer;
206 ✗ PetscViewerCreate(m_petscComm,&viewer);
207 ✗ PetscViewerASCIIOpen(m_petscComm,fileName.c_str(),&viewer);
208 ✗ PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
209 ✗ vector.view(viewer);
210 }
211
212
213 // Print solution of the system solve.
214 ✗ void CoefProblem::printSolution(int verbose, std::ostream& outstr) const {
215 IGNORE_UNUSED_OUTSTR;
216 ✗ if (verbose > 19) {
217 ✗ PetscPrintf(m_petscComm,"\n/================Solution of the assembly system Ax=b================/\n");
218 ✗ m_sol.view();
219 }
220 }
221
222 }
223