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 |