GCC Code Coverage Report


Directory: ./
File: Solver/linearProblemNSFracStepAdvDiff.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 147 160 91.9%
Branches: 86 145 59.3%

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/linearProblemNSFracStepAdvDiff.hpp"
21 #include "Core/felisceTransient.hpp"
22 #include "FiniteElement/elementVector.hpp"
23 #include "FiniteElement/elementMatrix.hpp"
24
25 namespace felisce {
26 8 LinearProblemNSFracStepAdvDiff::LinearProblemNSFracStepAdvDiff():
27 LinearProblem("Navier-Stokes Fractional Step: Adv-Diff"),
28 8 m_viscosity(0.),
29 8 m_density(0.),
30 8 m_explicitAdvection(0.),
31 8 m_buildTeporaryMatrix(false),
32 8 m_bdf(nullptr),
33 8 allocateSeqVec(false),
34 8 allocateSeqVecExt(false),
35
15/30
✓ 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.
✓ Branch 25 taken 8 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 8 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 8 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 8 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 8 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 8 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 8 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 8 times.
✗ Branch 47 not taken.
8 allocateSeqVecIncremental(false)
36 8 {}
37
38 16 LinearProblemNSFracStepAdvDiff::~LinearProblemNSFracStepAdvDiff() {
39
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
16 if(m_buildTeporaryMatrix)
40 8 m_matrix.destroy();
41 16 m_bdf = nullptr;
42
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
16 if (allocateSeqVec) {
43 16 m_seqBdfRHS.destroy();
44 }
45
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
16 if (allocateSeqVecExt) {
46 16 m_solExtrapol.destroy();
47 }
48 }
49
50 8 void LinearProblemNSFracStepAdvDiff::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) {
51
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 LinearProblem::initialize(mesh,comm, doUseSNES);
52 8 m_fstransient = fstransient;
53 8 std::vector<PhysicalVariable> listVariable;
54 8 std::vector<std::size_t> listNumComp;
55
56
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 listVariable.push_back(velocity);
57
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
8 listNumComp.push_back(this->dimension());
58
59
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 listVariable.push_back(pressure);
60
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 listNumComp.push_back(1);
61
62
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 > 0 ) {
63
2/10
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
4 FEL_CHECK(FelisceParam::instance().useALEformulation == 1,"Only standard ALE scheme implemented");
64
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
4 listVariable.push_back(displacement);
65
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
4 listNumComp.push_back(this->dimension());
66 }
67
68 //define unknown of the linear system.
69
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 m_listUnknown.push_back(velocity);
70
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 definePhysicalVariable(listVariable,listNumComp);
71
72
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 m_viscosity = FelisceParam::instance().viscosity;
73
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 m_density = FelisceParam::instance().density;
74
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 m_useSymmetricStress = FelisceParam::instance().useSymmetricStress;
75 8 }
76
77 164 void LinearProblemNSFracStepAdvDiff::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS) {
78 IGNORE_UNUSED_ELT_TYPE;
79 IGNORE_UNUSED_FLAG_MATRIX_RHS;
80
81 164 m_iVelocity = m_listVariable.getVariableIdList(velocity);
82 164 m_iPressure = m_listVariable.getVariableIdList(pressure);
83 164 m_feVel = m_listCurrentFiniteElement[m_iVelocity];
84 164 m_fePres = m_listCurrentFiniteElement[m_iPressure];
85 164 m_velocity = &m_listVariable[m_iVelocity];
86 164 m_pressure = &m_listVariable[m_iPressure];
87
88 /* if (m_explicitAdvection) {
89 m_iVelocityAdvection = m_listVariable.getVariableIdList(velocityAdvection);
90 m_feVelAdv = m_listCurrentFiniteElement[m_iVelocityAdvection];
91 m_velocityAdvection = &m_listVariable[m_iVelocityAdvection];
92 m_elemFieldAdv.initialize(DOF_FIELD,*m_feVelAdv,this->dimension());
93 }
94 */
95
96 164 m_elemFieldRHS.initialize(DOF_FIELD,*m_feVel,this->dimension());
97 164 m_elemFieldPres.initialize(DOF_FIELD,*m_fePres,1);
98 164 m_elemFieldPresDelta.initialize(DOF_FIELD,*m_fePres,1);
99
100 // Bdf
101 164 m_elemFieldRHSbdf.initialize(DOF_FIELD,*m_feVel,this->dimension());
102 164 m_elemFieldExt.initialize(DOF_FIELD,*m_feVel,this->dimension());
103
104 //=========
105 // ALE
106 //=========
107
2/2
✓ Branch 1 taken 120 times.
✓ Branch 2 taken 44 times.
164 if(FelisceParam::instance().useALEformulation > 0) {
108
109 120 m_iDisplacement = m_listVariable.getVariableIdList(displacement);
110 120 m_displacement = &m_listVariable[m_iDisplacement];
111
112 120 m_elemFieldVelMesh.initialize(DOF_FIELD,*m_feVel,this->dimension());
113
114 120 numDofExtensionProblem = m_externalDof[1]->numDof();
115
116 120 m_petscToGlobal1.resize(numDofExtensionProblem);
117 120 m_petscToGlobal2.resize(numDofExtensionProblem);
118 120 m_auxvec.resize(numDofExtensionProblem);
119
120
2/2
✓ Branch 0 taken 2547720 times.
✓ Branch 1 taken 120 times.
2547840 for (int i=0; i<numDofExtensionProblem; i++) {
121 2547720 m_petscToGlobal1[i]=i;
122 2547720 m_petscToGlobal2[i]=i;
123 }
124 120 AOApplicationToPetsc(m_externalAO[1],numDofExtensionProblem,m_petscToGlobal1.data());
125 120 AOApplicationToPetsc(m_ao,numDofExtensionProblem,m_petscToGlobal2.data());
126
127
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
120 if ( m_fstransient->iteration > 0) {
128
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 112 times.
120 if ( m_fstransient->iteration == 1) {
129 8 m_beta.duplicateFrom(m_solExtrapol);
130 }
131
132 //=========
133 // ALE
134 //=========
135 // Construction of m_beta = u^{n-1} - w^{n}
136 // Recall that externalVec(1) = - w^{n}
137 120 m_beta.copyFrom(m_solExtrapol);
138 120 externalVec(1).getValues(numDofExtensionProblem,m_petscToGlobal1.data(),m_auxvec.data());
139 120 m_beta.setValues(numDofExtensionProblem,m_petscToGlobal2.data(),m_auxvec.data(),ADD_VALUES);
140 }
141 }
142 164 }
143
144 1210894 void LinearProblemNSFracStepAdvDiff::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
145 IGNORE_UNUSED_FLAG_MATRIX_RHS;
146 IGNORE_UNUSED_ELEM_ID_POINT;
147
148 // note:
149 1210894 m_feVel->updateFirstDeriv(0, elemPoint);
150 1210894 m_fePres->updateFirstDeriv(0, elemPoint);
151
152 // in the NSFracStepExplicitModel, the unknowns are:
153 // velocityAdv[0], velocity[1], Pressure[2],
154 // otherwise: velocity[0], Pressure[1]/
155 //if (m_explicitAdvection) {
156 // m_feVelAdv->updateFirstDeriv(0, elemPoint);
157 //}
158
159 1210894 double coef = m_density/m_fstransient->timeStep;
160
161
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1210894 times.
1210894 if ( m_feVel->measure() < 0. ) {
162 int iglo;
163 ISLocalToGlobalMappingApply(this->m_mappingElem[m_currentMesh],1,&iel,&iglo);
164 FEL_ERROR("Error: Element of negative measure.");
165 }
166
167
168
5/6
✓ Branch 1 taken 5354 times.
✓ Branch 2 taken 1205540 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 5354 times.
✓ Branch 6 taken 1157354 times.
✓ Branch 7 taken 53540 times.
2416434 if ( (m_fstransient->iteration == 0 && FelisceParam::instance().useALEformulation == 0 ) ||
169
3/4
✓ Branch 1 taken 1205540 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1152000 times.
✓ Branch 5 taken 53540 times.
1205540 (m_fstransient->iteration > 0 && FelisceParam::instance().useALEformulation > 0 ) ) {
170
171
4/4
✓ Branch 0 taken 1152000 times.
✓ Branch 1 taken 5354 times.
✓ Branch 2 taken 576000 times.
✓ Branch 3 taken 576000 times.
1157354 if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix) ) {
172 //================================
173 // matrix
174 //================================
175
176
2/2
✓ Branch 0 taken 576000 times.
✓ Branch 1 taken 5354 times.
581354 if (m_useSymmetricStress) {
177 // \eta \int \epsilon u^{n+1} : \epsilon v
178 576000 m_elementMat[0]->eps_phi_i_eps_phi_j(2*m_viscosity,*m_feVel, 0, 0, this->dimension());
179 } else {
180 // \eta \int \nabla u^{n+1} \nabla v
181 5354 m_elementMat[0]->grad_phi_i_grad_phi_j(m_viscosity, *m_feVel, 0, 0, this->dimension());
182 }
183
184 // \dfrac{\rho}{\dt} u^{n+1} v
185 // 1/ no bdf :
186 // m_elementMat[0]->phi_i_phi_j(coef,*m_feVel,0,0,this->dimension());
187 // 2/ bdf :
188 581354 m_elementMat[0]->phi_i_phi_j(m_bdf->coeffDeriv0()*coef,*m_feVel,0,0,this->dimension());
189
190 //================================
191 // FE stabilization
192 //================================
193
194
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 581354 times.
581354 if (this->verbosity()>3) {
195 if ( m_velocity->finiteElementType() == m_pressure->finiteElementType())
196 std::cout << "The FE stabilization term has to be implemented in LinearProblemNSFracStepAdvDiff class" << std::endl;
197 }
198 }
199
200 }
201
202
2/2
✓ Branch 1 taken 1205540 times.
✓ Branch 2 taken 5354 times.
1210894 if ( m_fstransient->iteration > 0) {
203
204
4/4
✓ Branch 0 taken 1152000 times.
✓ Branch 1 taken 53540 times.
✓ Branch 2 taken 576000 times.
✓ Branch 3 taken 576000 times.
1205540 if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix) ) {
205
206
1/3
✗ Branch 1 not taken.
✓ Branch 2 taken 629540 times.
✗ Branch 3 not taken.
629540 switch(FelisceParam::instance().NSequationFlag) {
207 case 0:
208 //nothig to do: Stoke's model
209 break;
210
211 629540 case 1:
212 //================================
213 // matrix
214 //================================
215
216 //if (!m_explicitAdvection) {
217 // convective term
218
219
220 // 1/ no bdf
221 // m_elemFieldAdv.setValue(m_seqSol, *m_feVel, iel, m_iVelocity, m_ao, m_dof);
222 // m_elementMat[0]->u_grad_phi_j_phi_i(m_density,m_elemFieldAdv,*m_feVel,0,0,this->dimension());
223
224
225 // Temam's operator to stabilize the convective term
226 // --
227 // Convective term in energy balance :
228 // ... + \rho/2 \intGamma (u^n.n)|u^n+1|^2 - \rho/2 \intOmega (\nabla . u^n) |u^n+1|^2 + ...
229 // Add \rho/2 \intOmega (\nabla . u^n) u^n+1 v to delete the second term, which is consistent with the exact solution (divergence = 0)
230 // The flux of kinetic energy still remains...
231 // m_elementMat[0]->div_u_phi_j_phi_i(0.5,m_density,m_elemFieldExt,*m_feVel,0,0,this->dimension());
232 //}
233
234
235
2/2
✓ Branch 1 taken 576000 times.
✓ Branch 2 taken 53540 times.
629540 if(FelisceParam::instance().useALEformulation > 0) {
236 //=========
237 // ALE
238 //=========
239 // m_elemFieldAdv -> m_beta, where m_beta = u^{n-1} - w^{n} is defined in userElementInit()
240 576000 m_elemFieldExt.setValue(m_beta, *m_feVel, iel, m_iVelocity, m_ao, dof());
241 576000 m_elementMat[0]->u_grad_phi_j_phi_i(m_density,m_elemFieldExt,*m_feVel,0,0,this->dimension());
242
243 // m_elemFielVelMesh -> -w^{n}
244 576000 m_elemFieldVelMesh.setValue(externalVec(1), *m_feVel, iel, m_iDisplacement, m_externalAO[1], *m_externalDof[1]);
245 576000 m_elementMat[0]->div_u_phi_j_phi_i(m_density,m_elemFieldVelMesh,*m_feVel,0,0,this->dimension());
246 } else {
247 // 2/ bdf
248 // extrapolation of u^n in (u^n \cdot \grad{u^{n+1}})
249 53540 m_elemFieldExt.setValue(m_solExtrapol, *m_feVel, iel, m_iVelocity, m_ao, dof());
250 53540 m_elementMat[0]->u_grad_phi_j_phi_i(m_density,m_elemFieldExt,*m_feVel,0,0,this->dimension());
251 }
252
253 // SUPG stabilization
254 1259080 m_elementMat[0]->stab_supgAdvDiffCT(m_fstransient->timeStep, FelisceParam::instance().stabSUPG,
255 629540 FelisceParam::instance().stabdiv, m_density, m_viscosity, *m_feVel, m_elemFieldExt);
256 629540 break;
257 default:
258 FEL_ERROR("Error: No Convection Method found for NSFracStep solver");
259 }
260 }
261
262
4/4
✓ Branch 0 taken 1152000 times.
✓ Branch 1 taken 53540 times.
✓ Branch 2 taken 576000 times.
✓ Branch 3 taken 576000 times.
1205540 if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs) ) {
263 //================================
264 // RHS
265 //================================
266
267 //if (!m_explicitAdvection) {
268 // pressure term (non incremental version, but all the equation is written with the velocity of the prediction stage)
269 629540 m_elemFieldPres.setValue(externalVec(0), *m_fePres, iel, m_iPressure, m_externalAO[0], *m_externalDof[0]);
270
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 629540 times.
629540 assert(!m_elementVector.empty());
271 629540 m_elementVector[0]->grad_f_phi_i(-1,*m_feVel,*m_fePres,m_elemFieldPres,0);
272 // Incremental scheme
273
2/2
✓ Branch 1 taken 576000 times.
✓ Branch 2 taken 53540 times.
629540 if ( FelisceParam::instance().orderPressureExtrapolation > 0 ) { // && (m_fstransient->iteration > 1)
274 576000 m_elemFieldPresDelta.setValue(m_seqPressureDelta, *m_fePres, iel, m_iPressure, m_externalAO[0], *m_externalDof[0]);
275 576000 m_elementVector[0]->grad_f_phi_i(-1.,*m_feVel,*m_fePres,m_elemFieldPresDelta,0);
276 }
277
278 // time scheme bdf
279 // \frac{\rho}{\dt} \sum_{i=1}^k \alpha_i u_{n+1-i}
280 629540 m_elemFieldRHSbdf.setValue(m_seqBdfRHS, *m_feVel, iel, m_iVelocity, m_ao, dof());
281 629540 m_elementVector[0]->source(m_density,*m_feVel,m_elemFieldRHSbdf,0,this->dimension());
282 // 'm_density' only because the time step is integrated into the bdf term
283
284 // SUPG/PSPG stabilization
285
1/2
✓ Branch 1 taken 629540 times.
✗ Branch 2 not taken.
629540 if(FelisceParam::instance().NSequationFlag == 1) {
286
2/2
✓ Branch 1 taken 576000 times.
✓ Branch 2 taken 53540 times.
629540 if(FelisceParam::instance().useALEformulation > 0)
287 576000 m_elemFieldExt.setValue(m_beta, *m_feVel, iel, m_iVelocity, m_ao, dof());
288 else
289 53540 m_elemFieldExt.setValue(m_solExtrapol, *m_feVel, iel, m_iVelocity, m_ao, dof());
290 }
291 else
292 m_elemFieldExt.val.clear(); // zero convective velocity
293 1259080 m_elementVector[0]->stab_supgAdvDiffCT(m_fstransient->timeStep, FelisceParam::instance().stabSUPG, m_density, m_viscosity,
294 629540 *m_feVel, *m_fePres, m_elemFieldExt, m_elemFieldPres, m_elemFieldRHS);
295 }
296 }
297 1210894 }
298
299 void LinearProblemNSFracStepAdvDiff::computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
300 IGNORE_UNUSED_FLAG_MATRIX_RHS;
301 IGNORE_UNUSED_IEL;
302 IGNORE_UNUSED_ELEM_ID_POINT;
303 assert(!m_elementVectorBD.empty());
304
305 // General updates
306
307 for (std::size_t iFe = 0; iFe < m_listCurvilinearFiniteElement.size(); iFe++) {
308 m_listCurvilinearFiniteElement[iFe]->updateMeasNormal(0, elemPoint);
309 }
310
311 // Saving ptr to velocity curvilinear finite element
312 // CurvilinearFiniteElement* curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
313
314
315 }
316
317 8 void LinearProblemNSFracStepAdvDiff::initExtrapol(PetscVector& V_1) {
318
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
8 if (allocateSeqVecExt == false) {
319 8 m_solExtrapol.create(PETSC_COMM_SELF);
320 8 m_solExtrapol.setType(VECSEQ);
321 8 m_solExtrapol.setSizes(PETSC_DECIDE,m_numDof);
322 8 allocateSeqVecExt = true;
323 }
324 8 gatherVector(V_1, m_solExtrapol);
325 8 }
326
327 92 void LinearProblemNSFracStepAdvDiff::updateExtrapol(PetscVector& V_1) {
328 92 gatherVector(V_1, m_solExtrapol);
329 92 }
330
331 100 void LinearProblemNSFracStepAdvDiff::gatherVectorBeforeAssembleMatrixRHS() {
332
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 92 times.
100 if (allocateSeqVec == false) {
333 8 m_seqBdfRHS.create(PETSC_COMM_SELF);
334 8 m_seqBdfRHS.setType(VECSEQ);
335 8 m_seqBdfRHS.setSizes(PETSC_DECIDE,m_numDof);
336 8 allocateSeqVec = true;
337 }
338 100 gatherVector(m_bdf->vector(), m_seqBdfRHS);
339
340 // Incremental scheme
341
2/2
✓ Branch 1 taken 60 times.
✓ Branch 2 taken 40 times.
100 if (FelisceParam::instance().orderPressureExtrapolation > 0) {
342
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 56 times.
60 if (allocateSeqVecIncremental == false) {
343 4 m_seqPressureDelta.duplicateFrom(externalVec(0));
344 4 m_seqPressure0.duplicateFrom(externalVec(0));
345 4 m_seqPressure0.copyFrom(externalVec(0));
346 4 allocateSeqVecIncremental = true;
347 }
348
1/2
✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
60 waxpy(m_seqPressureDelta, -1., m_seqPressure0, externalVec(0));
349 60 m_seqPressure0.copyFrom(externalVec(0));
350 }
351
352 100 }
353
354 4 void LinearProblemNSFracStepAdvDiff::copyMatrixRHS() {
355 4 m_matrix.duplicateFrom(matrix(0),MAT_COPY_VALUES);
356 4 m_matrix.assembly(MAT_FINAL_ASSEMBLY);
357 4 m_buildTeporaryMatrix = true;
358 4 }
359
360 40 void LinearProblemNSFracStepAdvDiff::addMatrixRHS() {
361 40 matrix(0).axpy(1,m_matrix,SAME_NONZERO_PATTERN);
362 40 }
363 }
364