GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemNSFracStepProj.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 144 155 92.9%
Branches: 105 216 48.6%
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-F. Gerbeau
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/linearProblemNSFracStepProj.hpp"
21 #include "FiniteElement/elementVector.hpp"
22 #include "FiniteElement/elementMatrix.hpp"
23 #include "DegreeOfFreedom/boundaryCondition.hpp"
24
25 namespace felisce {
26 8 LinearProblemNSFracStepProj::LinearProblemNSFracStepProj():
27 LinearProblem("Navier-Stokes Fractional Step: Projection"),
28 8 m_viscosity(0.),
29 8 m_density(0.),
30 8 m_lumpedModelBC(nullptr),
31
7/14
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 8 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 8 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 8 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 8 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 8 times.
✗ Branch 23 not taken.
 8 m_buildTeporaryMatrix(false)
32 8 {}
33
34 16 LinearProblemNSFracStepProj::~LinearProblemNSFracStepProj() {
35
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 16 if(m_buildTeporaryMatrix)
36 8 m_matrix.destroy();
37 }
38
39 8 void LinearProblemNSFracStepProj::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) {
40
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 LinearProblem::initialize(mesh,comm, doUseSNES);
41 8 m_fstransient = fstransient;
42 8 std::vector<PhysicalVariable> listVariable;
43 8 std::vector<std::size_t> listNumComp;
44
45
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_explicitAdvection = FelisceParam::instance().explicitAdvection;
46
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 if (m_explicitAdvection) {
47 ✗ listVariable.push_back(velocityAdvection);
48 ✗ listNumComp.push_back(this->dimension());
49 }
50
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 listVariable.push_back(velocity);
51
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 listNumComp.push_back(this->dimension());
52
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 listVariable.push_back(pressure);
53
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 listNumComp.push_back(1);
54
55
3/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 4 times.
 8 if(FelisceParam::instance().useALEformulation) {
56
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listVariable.push_back(displacement);
57
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 listNumComp.push_back(this->dimension());
58 }
59
60 //define unknown of the linear system.
61
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_listUnknown.push_back(pressure);
62
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 definePhysicalVariable(listVariable,listNumComp);
63
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_viscosity = FelisceParam::instance().viscosity;
64
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_density = FelisceParam::instance().density;
65 8 }
66
67 4 void LinearProblemNSFracStepProj::initializeLumpedModelBC(LumpedModelBC* lumpedModelBC) {
68 4 m_lumpedModelBC = lumpedModelBC;
69 4 }
70
71 104 void LinearProblemNSFracStepProj::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS) {
72 IGNORE_UNUSED_ELT_TYPE;
73 IGNORE_UNUSED_FLAG_MATRIX_RHS;
74 104 m_iVelocity = m_listVariable.getVariableIdList(velocity);
75 104 m_iPressure = m_listVariable.getVariableIdList(pressure);
76 104 m_feVel = m_listCurrentFiniteElement[m_iVelocity];
77 104 m_fePres = m_listCurrentFiniteElement[m_iPressure];
78 104 m_velocity = &m_listVariable[m_iVelocity];
79 104 m_pressure = &m_listVariable[m_iPressure];
80
81 104 m_elemFieldAdv.initialize(DOF_FIELD,*m_feVel,this->dimension());
82 104 m_elemFieldRHS.initialize(DOF_FIELD,*m_feVel,this->dimension());
83 104 m_elemFieldDOF.initialize(DOF_FIELD,*m_feVel,this->dimension());
84 104 m_elemFieldPress.initialize(DOF_FIELD,*m_fePres,1);
85 104 }
86
87 634894 void LinearProblemNSFracStepProj::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
88 IGNORE_UNUSED_FLAG_MATRIX_RHS;
89 IGNORE_UNUSED_IEL;
90 IGNORE_UNUSED_ELEM_ID_POINT;
91
92 634894 m_feVel->updateFirstDeriv(0, elemPoint);
93 634894 m_fePres->updateFirstDeriv(0, elemPoint);
94 634894 double coef = m_fstransient->timeStep/m_density;
95
96
97
98
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 634894 times.
 634894 if ( m_feVel->measure() < 0. ) {
99 int iglo;
100 ✗ ISLocalToGlobalMappingApply(this->m_mappingElem[m_currentMesh],1,&iel,&iglo);
101 ✗ std::cout << " iglo = " << iglo << std::endl;
102 ✗ FEL_ERROR("Error: Element of negative measure.");
103 }
104
105
106 //================================
107 // matrix
108 //================================
109
5/6
✓ Branch 1 taken 5354 times.
✓ Branch 2 taken 629540 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 5354 times.
✓ Branch 6 taken 581354 times.
✓ Branch 7 taken 53540 times.
 1264434 if ( (m_fstransient->iteration == 0 && FelisceParam::instance().useALEformulation == 0 ) ||
110
3/4
✓ Branch 1 taken 629540 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 576000 times.
✓ Branch 5 taken 53540 times.
 629540 (m_fstransient->iteration > 0 && FelisceParam::instance().useALEformulation > 0 ) ) {
111
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 581354 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 581354 if( flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_matrix )
112 // - \int \frac{dt}{rhof} \nabla p^{n} \nabla q
113 581354 m_elementMat[0]->grad_phi_i_grad_phi_j(coef,*m_fePres,0,0,1);
114 }
115
116
2/2
✓ Branch 1 taken 629540 times.
✓ Branch 2 taken 5354 times.
 634894 if ( m_fstransient->iteration > 0 ) {
117 //SUPG/PSPG Stabilization
118
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 629540 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 629540 if( flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_matrix ) {
119 629540 m_elemFieldAdv.setValue(externalVec(1), *m_feVel, iel, m_iVelocity, m_externalAO[0], *m_externalDof[0]);
120 629540 m_elementMat[0]->stab_supgProjCT(m_fstransient->timeStep, FelisceParam::instance().stabSUPG, m_density, m_viscosity, *m_fePres, m_elemFieldAdv);
121 }
122 }
123
124 //================================
125 // RHS
126 //================================
127
2/2
✓ Branch 1 taken 629540 times.
✓ Branch 2 taken 5354 times.
 634894 if (m_fstransient->iteration > 0) {
128
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 629540 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 629540 if( flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_rhs ) {
129 // - \int div \tilde{u}^n
130 629540 m_elemFieldDOF.setValue(externalVec(0), *m_feVel, iel, m_iVelocity, m_externalAO[0], *m_externalDof[0]);
131
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 629540 times.
 629540 assert(!m_elementVector.empty());
132 629540 m_elementVector[0]->div_f_phi_i(-1.,*m_feVel,*m_fePres,m_elemFieldDOF,0);
133
134 // SUPG/PSPG Stabilization
135 629540 m_elemFieldAdv.setValue(externalVec(1), *m_feVel, iel, m_iVelocity, m_externalAO[0], *m_externalDof[0]);
136 629540 m_elemFieldPress.setValue(this->sequentialSolution(), *m_fePres, iel, m_iPressure, m_ao, dof());
137 629540 m_elementVector[0]->stab_supgProjCT(m_fstransient->timeStep, FelisceParam::instance().stabSUPG, m_density, m_viscosity, *m_fePres, m_elemFieldAdv, m_elemFieldDOF, m_elemFieldPress, m_elemFieldRHS);
138 }
139 }
140 634894 }
141
142
143 108 void LinearProblemNSFracStepProj::finalizeEssBCTransientDerivedProblem() {
144
2/2
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 64 times.
 108 if (m_lumpedModelBC) {
145
1/2
✓ Branch 2 taken 44 times.
✗ Branch 3 not taken.
 44 if(FelisceParam::instance().model == "NSFracStep") {
146 BoundaryCondition* bc;
147 44 std::vector<double> PlumpedModelBC;
148
2/2
✓ Branch 1 taken 88 times.
✓ Branch 2 taken 44 times.
 132 for(std::size_t iLumpedModelBC=0; iLumpedModelBC<m_lumpedModelBC->size(); iLumpedModelBC++) {
149
2/4
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
 88 if(FelisceParam::instance().lumpedModelBCAlgo[iLumpedModelBC] == 1) {
150 // get the boundary condition we want to std::set the value
151 88 bc = getBoundaryConditionList().EssentialLumpedModelBC(iLumpedModelBC);
152
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 88 times.
 88 if ( m_linearizedFlag )
153 ✗ PlumpedModelBC.push_back(0.0);
154 else {
155
2/4
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 88 times.
 88 if ( FelisceParam::instance().orderPressureExtrapolation > 0 )
156 ✗ PlumpedModelBC.push_back(this->m_lumpedModelBC->Pp(iLumpedModelBC)-this->m_lumpedModelBC->PpOld(iLumpedModelBC));
157 else
158
1/2
✓ Branch 2 taken 88 times.
✗ Branch 3 not taken.
 88 PlumpedModelBC.push_back(this->m_lumpedModelBC->Pp(iLumpedModelBC));
159 }
160
161
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 88 times.
 88 if (this->verbosity()>2) {
162 ✗ std::cout << "LumpedModelBC :" << std::endl;
163 ✗ std::cout << "Label = " << m_currentLabel << std::endl;
164 ✗ std::cout << "Applied PlumpedModelBC = " << PlumpedModelBC[0] << std::endl;
165 }
166
167 // std::set boundary condition
168
1/2
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
 88 bc->setValue(PlumpedModelBC);
169 88 PlumpedModelBC.clear();
170 }
171 }
172 44 }
173 }
174 108 }
175
176 100 void LinearProblemNSFracStepProj::applyEssentialBoundaryConditionDerivedProblem(int rank, FlagMatrixRHS flagMatrixRHS) {
177
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 60 times.
 100 if (m_lumpedModelBC) {
178
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 if(FelisceParam::instance().model == "NSFracStep") {
179 BoundaryCondition* BC;
180 felInt ia;
181 40 std::unordered_set<felInt> idDofBC;
182 40 std::unordered_map<int, double> idBCAndValue;
183 40 std::unordered_set<felInt> idDofOfAllEssLumpedModelBC;
184
185 // non symmetric pseudo-elimination, cf. LinearProblem::applyEssentialBoundaryConditionHelper<DirichletApplicationOptions::nonSymmetricPseudoElimination>
186
2/2
✓ Branch 1 taken 80 times.
✓ Branch 2 taken 40 times.
 120 for (std::size_t i=0; i < m_boundaryConditionList.numEssentialLumpedModelBoundaryCondition(); i++) {
187
2/4
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 80 times.
✗ Branch 5 not taken.
 80 if(FelisceParam::instance().lumpedModelBCAlgo[i] == 1) {
188 80 BC = m_boundaryConditionList.EssentialLumpedModelBC(i);
189 80 idBCAndValue.clear();
190 80 idDofBC.clear();
191
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 getListDofOfBC(*BC, idDofBC, idBCAndValue);
192
2/2
✓ Branch 4 taken 400 times.
✓ Branch 5 taken 80 times.
 480 for(auto itDofBC = idDofBC.begin(); itDofBC != idDofBC.end(); itDofBC++) {
193
1/2
✓ Branch 2 taken 400 times.
✗ Branch 3 not taken.
 400 idDofOfAllEssLumpedModelBC.insert(*itDofBC);
194 400 ia = *itDofBC;
195
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 400 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 400 if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs))
196
2/4
✓ Branch 3 taken 400 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 400 times.
✗ Branch 7 not taken.
 400 vector().setValue( ia, idBCAndValue[*itDofBC], INSERT_VALUES);
197 }
198 }
199 }
200
201 // get the number of dof in each process for all EssentialLumpedModel BC
202 int numProc;
203
2/4
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40 times.
✗ Branch 5 not taken.
 40 MPI_Comm_size(MpiInfo::petscComm(), &numProc);
204
2/4
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 40 times.
✗ Branch 5 not taken.
 40 int *sizeTmp = new int [numProc];
205
2/4
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 40 times.
✗ Branch 5 not taken.
 40 int *sizeDofBC = new int [numProc];
206
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 40 times.
 200 for (int rankProc = 0; rankProc < numProc; rankProc++) {
207 160 sizeTmp[rankProc] = 0;
208 }
209 40 sizeTmp[rank] = idDofOfAllEssLumpedModelBC.size();
210
2/4
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40 times.
✗ Branch 5 not taken.
 40 MPI_Allreduce(&sizeTmp[0], &sizeDofBC[0], numProc, MPI_INT, MPI_SUM, MpiInfo::petscComm());
211
212 40 int startDof = 0;
213
2/2
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 40 times.
 100 for (int rankProc = 0; rankProc < rank; rankProc++) {
214 60 startDof += sizeDofBC[rankProc];
215 }
216 // get the number of dof in all processes = sum of sizeDofBC[] (some dof are counted several times)
217 40 int numDofBCTotal = 0;
218
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 40 times.
 200 for (int rankProc = 0; rankProc < numProc; rankProc++) {
219 160 numDofBCTotal += sizeDofBC[rankProc];
220 }
221
222 // put all the id dof in an array (some dof are counted several times)
223
2/4
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 40 times.
✗ Branch 5 not taken.
 40 int *dofBCTmp = new int [numDofBCTotal];
224
2/4
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 40 times.
✗ Branch 5 not taken.
 40 int *dofBCArray = new int [numDofBCTotal];
225
2/2
✓ Branch 0 taken 1600 times.
✓ Branch 1 taken 40 times.
 1640 for (int iDof = 0; iDof < numDofBCTotal; iDof++) {
226 1600 dofBCTmp[iDof] = 0;
227 }
228 40 int iDof = 0;
229
2/2
✓ Branch 3 taken 400 times.
✓ Branch 4 taken 40 times.
 440 for(auto itDofBC = idDofOfAllEssLumpedModelBC.begin(); itDofBC != idDofOfAllEssLumpedModelBC.end(); ++itDofBC) {
230 400 dofBCTmp[startDof + iDof] = *itDofBC;
231 400 iDof++;
232 }
233
2/4
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40 times.
✗ Branch 5 not taken.
 40 MPI_Allreduce(&dofBCTmp[0], &dofBCArray[0], numDofBCTotal, MPI_INT, MPI_SUM, MpiInfo::petscComm());
234
235 // rearrange the id dof in the increasing order (each dof is counted only once)
236 40 idDofOfAllEssLumpedModelBC.clear();
237
2/2
✓ Branch 0 taken 1600 times.
✓ Branch 1 taken 40 times.
 1640 for (iDof = 0; iDof < numDofBCTotal; iDof++) {
238
1/2
✓ Branch 1 taken 1600 times.
✗ Branch 2 not taken.
 1600 idDofOfAllEssLumpedModelBC.insert(dofBCArray[iDof]);
239 }
240 40 numDofBCTotal = idDofOfAllEssLumpedModelBC.size();
241 40 iDof = 0;
242
2/2
✓ Branch 3 taken 1600 times.
✓ Branch 4 taken 40 times.
 1640 for(auto itDofBC = idDofOfAllEssLumpedModelBC.begin(); itDofBC != idDofOfAllEssLumpedModelBC.end(); ++itDofBC) {
243 1600 dofBCArray[iDof] = *itDofBC;
244 1600 iDof++;
245 }
246
247
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 40 if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix)) {
248
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 matrix(0).setOption(MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
249
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 matrix(0).zeroRows(numDofBCTotal, &dofBCArray[0], 1.);
250 }
251
252
1/2
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
 40 delete [] sizeTmp;
253
1/2
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
 40 delete [] sizeDofBC;
254
1/2
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
 40 delete [] dofBCTmp;
255
1/2
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
 40 delete [] dofBCArray;
256 40 }
257 }
258 100 }
259
260 4 void LinearProblemNSFracStepProj::copyMatrixRHS() {
261 4 m_matrix.duplicateFrom(matrix(0),MAT_COPY_VALUES);
262 4 m_matrix.assembly(MAT_FINAL_ASSEMBLY);
263 4 m_buildTeporaryMatrix = true;
264 4 }
265
266 40 void LinearProblemNSFracStepProj::addMatrixRHS() {
267 40 matrix(0).axpy(1,m_matrix,SAME_NONZERO_PATTERN);
268 40 }
269
270 //Embedded interface
271 ✗ void LinearProblemNSFracStepProj::userAssembleMatrixRHSEmbeddedInterface() {}
272
273 }
274