GCC Code Coverage Report

Directory: ./
File: Model/NSFracStepModel.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 170 364 46.7%
Branches: 124 352 35.2%
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:
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Model/NSFracStepModel.hpp"
21 #include "Solver/linearProblemHarmonicExtension.hpp"
22 #include "Solver/linearProblemNSFracStepAdvDiff.hpp"
23 #include "Solver/linearProblemNSFracStepProj.hpp"
24
25 namespace felisce {
26
17/34
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 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 29 taken 8 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 8 times.
✗ Branch 33 not taken.
✓ Branch 35 taken 8 times.
✗ Branch 36 not taken.
✓ Branch 38 taken 8 times.
✗ Branch 39 not taken.
✓ Branch 41 taken 8 times.
✗ Branch 42 not taken.
✓ Branch 44 taken 8 times.
✗ Branch 45 not taken.
✓ Branch 47 taken 8 times.
✗ Branch 48 not taken.
✓ Branch 50 taken 8 times.
✗ Branch 51 not taken.
✓ Branch 53 taken 8 times.
✗ Branch 54 not taken.
✓ Branch 56 taken 8 times.
✗ Branch 57 not taken.
 8 NSFracStepModel::NSFracStepModel():Model() {
27 8 m_lumpedModelBC = nullptr;
28
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_name = "NSFracStep";
29 8 m_matrixWithoutBCcreated = false;
30 8 }
31
32
33 32 NSFracStepModel::~NSFracStepModel() {
34
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 16 if(m_lumpedModelBC)
35
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 8 delete m_lumpedModelBC;
36 16 m_solExtrapolate.destroy();
37 32 }
38
39 //***************************** Initialization ************************************
40
41 8 void NSFracStepModel::initializeDerivedModel() {
42
43
2/2
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 4 times.
 8 if( FelisceParam::instance().lumpedModelBCLabel.size() ) {
44
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 m_lumpedModelBC = new LumpedModelBC(m_fstransient);
45 4 static_cast<LinearProblemNSFracStepProj*>(m_linearProblem[1])->initializeLumpedModelBC(m_lumpedModelBC);
46 }
47
48
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 if(FelisceParam::instance().useALEformulation ) {
49 // Resizing the auxiliary vectors needed in updateMesh()
50 4 m_numDofExtensionProblem = m_linearProblem[2]->numDof();
51
1/2
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
 4 m_dispMeshVector.resize( m_linearProblem[2]->supportDofUnknown(0).listNode().size() * m_linearProblem[2]->meshLocal()->numCoor() );
52 4 m_petscToGlobal.resize(m_numDofExtensionProblem);
53
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( FelisceParam::instance().readPreStressedDisplFromFile )
54 ✗ m_preStressedDispMeshVector.resize(m_numDofExtensionProblem);
55 // From PETSc order to global
56
2/2
✓ Branch 0 taken 84924 times.
✓ Branch 1 taken 4 times.
 84928 for (int ipg=0; ipg< m_numDofExtensionProblem; ipg++)
57 84924 m_petscToGlobal[ipg] = ipg;
58
59 4 AOApplicationToPetsc(m_linearProblem[2]->ao(),m_numDofExtensionProblem,m_petscToGlobal.data());
60 4 static_cast<LinearProblemNSFracStepAdvDiff*>(m_linearProblem[0])->pressureFWI() = 0;
61
62
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
 4 if ( FelisceParam::instance().readDisplFromFile || FelisceParam::instance().readPreStressedDisplFromFile ) {
63 // defined in linearProblemHarmonicExtension
64 ✗ LinearProblemHarmonicExtension* lpHarmExt = static_cast<LinearProblemHarmonicExtension*>(m_linearProblem[2]);
65 ✗ lpHarmExt->readMatchFile_ALE(m_linearProblem[2]->listVariable(), FelisceParam::instance().MoveMeshMatchFile);
66 ✗ lpHarmExt->identifyIdDofToMatch_ALE(m_linearProblem[2]->dof());
67
68 ✗ if (FelisceParam::instance().hasInitialCondition || FelisceParam::instance().readPreStressedDisplFromFile ) {
69 ✗ PetscPrintf(MpiInfo::petscComm(),"\nCycle 1, Read data displacement at time t = 0.00000\n");
70 ✗ lpHarmExt->readDataDisplacement(m_ios, 0.0);
71 }
72 }
73 }
74
75
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
 8 if (FelisceParam::instance().orderBdfFS > 2)
76 ✗ FEL_ERROR("BDF not yet implemented for order greater than 2 with fractional step model.");
77
78 8 m_bdfFS.defineOrder(FelisceParam::instance().orderBdfFS);
79 8 m_linearProblem[0]->initializeTimeScheme(&m_bdfFS);
80
81
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
 8 if ( FelisceParam::instance().fsiCoupling) {
82 ✗ m_rhsWihoutBC.duplicateFrom(m_linearProblem[0]->vector());
83 ✗ m_rhsPressWihoutBC.duplicateFrom(m_linearProblem[1]->vector());
84 ✗ m_viscousResidual.duplicateFrom(m_linearProblem[0]->vector());
85 ✗ m_residual.duplicateFrom(m_linearProblem[0]->vector());
86 ✗ m_lumpedMassRN.duplicateFrom(m_linearProblem[0]->vector());
87 ✗ m_seqResidual.duplicateFrom(m_linearProblem[0]->sequentialSolution());
88 ✗ m_seqViscousResidual.duplicateFrom(m_linearProblem[0]->sequentialSolution());
89 ✗ m_seqPressLoad.duplicateFrom(m_linearProblem[0]->sequentialSolution());
90 }
91
92
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 if (FelisceParam::instance().orderPressureExtrapolation > 0) {
93 4 m_prevPress.duplicateFrom(m_linearProblem[1]->solution());
94 4 m_prevPress.zeroEntries();
95 }
96
97 8 }
98
99
100 8 void NSFracStepModel::setExternalVec() {
101 // std::set external vectors for linear problems
102 // (to be used in assembling routines)
103
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
 8 if (FelisceParam::instance().explicitAdvection) {
104 // this has been moved into NSFracStepExplicitModel
105 ✗ FEL_ERROR("NSFracStepModel should be used without explicitAdvection");
106
107 } else {
108 // standard fracstep
109 8 m_linearProblem[0]->pushBackExternalVec(m_linearProblem[1]->sequentialSolution()); // externalVec(0) in m_linearProblem[0]
110 8 m_linearProblem[1]->pushBackExternalVec(m_linearProblem[0]->sequentialSolution()); // externalVec(0) in m_linearProblem[1]
111
112 8 m_linearProblem[0]->pushBackExternalAO(m_linearProblem[1]->ao());
113 8 m_linearProblem[1]->pushBackExternalAO(m_linearProblem[0]->ao());
114
115 8 m_linearProblem[0]->pushBackExternalDof(m_linearProblem[1]->dof());
116 8 m_linearProblem[1]->pushBackExternalDof(m_linearProblem[0]->dof());
117
118 // Initialization of the std::vector containing the convective velocity
119 8 m_beta.duplicateFrom(m_linearProblem[0]->sequentialSolution());
120 8 m_beta.zeroEntries();
121
122 8 m_linearProblem[1]->pushBackExternalVec(m_beta); // externalVec(1) in m_linearProblem[1]
123
124
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 if(FelisceParam::instance().useALEformulation > 0) {
125
126 // fsi challenge
127
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( FelisceParam::instance().useElasticExtension ) {
128 ✗ m_dispMesho_pp.duplicateFrom(m_linearProblem[2]->solution());
129 ✗ m_dispMesho_pp.zeroEntries();
130 }
131
132 //1) Initialization of the std::vector containing the displacement of the mesh d^n
133 4 m_dispMesh.duplicateFrom(m_linearProblem[2]->sequentialSolution());
134 4 m_dispMesh.zeroEntries();
135
136 //2) Initialization of the std::vector containing the velocity of the mesh w^n
137 4 m_velMesh.duplicateFrom(m_linearProblem[2]->sequentialSolution());
138 4 m_velMesh.zeroEntries();
139
140 // Passing the ao and dof
141 4 m_linearProblem[0]->pushBackExternalAO(m_linearProblem[2]->ao()); //m_externalAO[1] in m_linearProblem[0]
142 4 m_linearProblem[0]->pushBackExternalDof(m_linearProblem[2]->dof()); //m_externalDof[1] in m_linearProblem[0]
143 // Passing the velocity of the mesh to the AdvDiff problem
144 4 m_linearProblem[0]->pushBackExternalVec(m_velMesh); //_externalVec[1] in m_linearProblem[0]
145 // Resizing the auxiliary vectors needed in updateMesh()
146 4 m_tmpDisp.resize(m_numDofExtensionProblem);
147
148
149
150 4 m_displTimeStep = 0.;
151 }
152
153 }
154
155 8 }
156
157
158 //************************** Pre & Post assembleMatrixRHS *******************************
159
160 20 void NSFracStepModel::preAssemblingMatrixRHS(std::size_t iProblem) {
161
162
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 12 times.
 20 if (iProblem == 0) {
163 8 m_solExtrapolate.duplicateFrom(m_linearProblem[0]->vector());
164 8 m_solExtrapolate.zeroEntries();
165
166 8 m_linearProblem[0]->solution().zeroEntries();
167 }
168
169
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 12 times.
 20 if(FelisceParam::instance().useALEformulation == 0) {
170 //First assembling loop in iteration 0 to build static matrix.
171 8 m_linearProblem[iProblem]->assembleMatrixRHS(MpiInfo::rankProc());
172 // save static matrix in matrix _A.
173 8 m_linearProblem[iProblem]->copyMatrixRHS();
174 } else {
175
5/6
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 8 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 8 times.
 12 if ( ( iProblem == 2 ) && ( !FelisceParam::instance().useElasticExtension ) ) {
176 // First assembling loop in iteration 0 to build static matrix of the extension problem
177 4 m_linearProblem[2]->assembleMatrixRHS(MpiInfo::rankProc());
178 // Save the static part of the extension problem in the matrix _A
179 4 m_linearProblem[2]->copyMatrixRHS();
180 }
181 }
182 20 }
183
184 260 void NSFracStepModel::postAssemblingMatrixRHS(std::size_t iProblem) {
185 // _Matrix = _A + _Matrix (add static matrix to the dynamic matrix to build
186 // complete matrix of the system
187
2/2
✓ Branch 1 taken 80 times.
✓ Branch 2 taken 180 times.
 260 if(FelisceParam::instance().useALEformulation == 0) {
188 80 m_linearProblem[iProblem]->addMatrixRHS();
189 } else {
190
2/2
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 120 times.
 180 if (iProblem == 2) {
191 60 m_linearProblem[2]->addMatrixRHS();
192 }
193 }
194 260 }
195
196
197 ✗ void NSFracStepModel::solveProjectionStep(bool linearFlag) {
198
199 ✗ m_linearProblem[1]->linearizedFlag() = linearFlag;
200
201 ✗ if ( linearFlag )
202 ✗ m_linearProblem[1]->vector().zeroEntries();
203 else
204 ✗ m_linearProblem[1]->vector().copyFrom(m_rhsPressWihoutBC);
205
206 ✗ if ( linearFlag && FelisceParam::verbose() )
207 ✗ PetscPrintf(MpiInfo::petscComm()," -- Projection-step linearized solver \n");
208 ✗ else if ( FelisceParam::verbose() )
209 ✗ PetscPrintf(MpiInfo::petscComm()," -- Projection-step solver \n");
210
211 ✗ if (FelisceParam::verbose())
212 ✗ PetscPrintf(MpiInfo::petscComm(),"Applying BC \n");
213
214 //Apply essential transient boundary conditions.
215 ✗ m_linearProblem[1]->finalizeEssBCTransient();
216
217 ✗ m_linearProblem[1]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs,
218 0, true, ApplyNaturalBoundaryConditions::yes);
219
220 ✗ PetscPrintf(MpiInfo::petscComm(),"Solving Projection-step (same matrix) \n");
221
222 //Solve linear system.
223 ✗ m_linearProblem[1]->solveWithSameMatrix();
224
225 ✗ if( (!linearFlag ) && (FelisceParam::instance().orderPressureExtrapolation > 0) ) // && (m_fstransient->iteration > 1)
226 ✗ m_linearProblem[1]->solution().axpy(1, m_prevPress);
227
228 ✗ m_linearProblem[1]->gatherSolution();
229
230 ✗ if (FelisceParam::verbose())
231 ✗ PetscPrintf(MpiInfo::petscComm(),"Evaluating fluid load \n");
232
233
234 // Pressure load calculation in FWI fashion
235 ✗ static_cast<LinearProblemNSFracStepAdvDiff*>(m_linearProblem[0])->pressureFWI() = 1;
236 ✗ m_linearProblem[0]->vector().set(0.);
237 ✗ m_linearProblem[0]->assembleMatrixRHSNaturalBoundaryCondition(MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, 0);
238 ✗ static_cast<LinearProblemNSFracStepAdvDiff*>(m_linearProblem[0])->pressureFWI() = 0;
239
240 ✗ if ( linearFlag )
241 ✗ m_residual.copyFrom(m_linearProblem[0]->vector());
242 else
243 ✗ waxpy(m_residual,1., m_linearProblem[0]->vector(),m_viscousResidual);
244
245 ✗ m_linearProblem[0]->gatherVector(m_residual, m_seqResidual);
246 ✗ m_linearProblem[0]->linearizedFlag() = false;
247
248 }
249
250 100 void NSFracStepModel::forward() {
251
252
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 60 times.
 160 if ( FelisceParam::instance().useALEformulation && FelisceParam::instance().initALEDispByFile
253
3/6
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 40 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 100 times.
 160 && (m_fstransient->iteration == 0 ) ) {
254
255 // init diplacement from file
256 ✗ std::vector<double> vector;
257 ✗ std::vector<double> vv;
258 ✗ std::vector<int> ix;
259 ✗ for (int ilinPb =1; ilinPb<=2; ++ilinPb)
260 {
261 int numDof;
262 ✗ if( FelisceParam::instance().FusionDof && (ilinPb == 2) )
263 ✗ numDof = m_linearProblem[ilinPb]->supportDofUnknown(0).listNode().size() * m_linearProblem[ilinPb]->meshLocal()->numCoor();
264 ✗ else if ( FelisceParam::instance().FusionDof && (ilinPb == 1) )
265 ✗ numDof = m_linearProblem[ilinPb]->supportDofUnknown(0).listNode().size();
266 else
267 ✗ numDof = m_linearProblem[ilinPb]->numDof();
268
269 ✗ vector.resize( numDof, 0.);
270
271 ✗ m_ios[0]->readVariable(ilinPb, 0.0, vector.data(), numDof);
272
273 ✗ if( FelisceParam::instance().FusionDof && (ilinPb == 1) ) {
274 ✗ int numNodes = m_linearProblem[ilinPb]->supportDofUnknown(0).listNode().size();
275 int idof;
276 ✗ std::vector<double> aux(vector);
277 ✗ for (int ip = 0; ip < numNodes; ip++) {
278 ✗ idof = m_linearProblem[1]->supportDofUnknown(0).listPerm()[ip];
279 ✗ vector[idof] = aux[ip];
280 }
281 }
282
283 ✗ if( FelisceParam::instance().FusionDof && (ilinPb == 2) ) {
284 ✗ int numCompDisp = m_linearProblem[2]->dimension();
285 ✗ int numCoorMesh = m_linearProblem[2]->meshLocal()->numCoor();
286 ✗ int numNodes = m_linearProblem[ilinPb]->supportDofUnknown(0).listNode().size();
287 int idof;
288 ✗ std::vector<double> aux(vector);
289 ✗ for (int ip = 0; ip < numNodes ; ip++) {
290 ✗ idof = m_linearProblem[2]->supportDofUnknown(0).listPerm()[ip];
291 ✗ for (felInt ic=0; ic < numCompDisp; ++ic)
292 ✗ vector[idof*numCompDisp+ic] = aux[ip*numCoorMesh+ic];
293 }
294 }
295
296 int low, high;
297 ✗ VecGetOwnershipRange(m_linearProblem[ilinPb]->solution().toPetsc(),&low,&high);
298
299 ✗ int localSize = high-low;
300
301 ✗ ix.resize(localSize);
302 ✗ int j=0;
303 ✗ for (int i=low; i<high;++i)
304 {
305 ✗ ix[j] = i;
306 ✗ ++j;
307 }
308 ✗ vv.resize(localSize);
309
310 ✗ AOPetscToApplication(m_linearProblem[ilinPb]->ao(),localSize,ix.data());
311 ✗ for (int i=0; i<localSize;++i)
312 ✗ vv[i] = vector[ix[i]];
313 ✗ AOApplicationToPetsc(m_linearProblem[ilinPb]->ao(),localSize,ix.data());
314
315 ✗ m_linearProblem[ilinPb]->solution().setValues(localSize,ix.data(),vv.data(),INSERT_VALUES);
316 ✗ m_linearProblem[ilinPb]->solution().assembly();
317 ✗ m_linearProblem[ilinPb]->gatherSolution();
318
319 }
320 }
321
322 // TO DO: pressure initialization challenge FSI to me removed
323 // if ( FelisceParam::instance().fsiCoupling && ( (FelisceParam::instance().testCase == 20) ||| (FelisceParam::instance().testCase == 23) )
324 // && FelisceParam::instance().useALEformulation && (m_fstransient->iteration == 0 ) ) {
325 // PetscInt localSize = m_linearProblem[1]->numDofLocalPerUnknown(m_linearProblem[1]->listUnknown()[0]);
326 // std::vector<felInt> idLocal(localSize);
327 // std::vector<felInt> idGlobalP(localSize);
328 // std::vector<felInt> idGlobalA(localSize);
329 // std::vector<double> value(localSize);
330 // for (felInt i = 0; i < localSize; ++i) {
331 // idLocal[i] = i;
332 // value[i]=0.;
333 // }
334 // ISLocalToGlobalMappingApply(m_linearProblem[1]->mappingDofLocalToDofGlobal(m_linearProblem[1]->listUnknown()[0]),localSize,idLocal.data(),idGlobalP.data());
335 // for (felInt i = 0; i < localSize; ++i)
336 // idGlobalA[i] = idGlobalP[i];
337
338 // AOPetscToApplication(m_linearProblem[1]->ao(),localSize,idGlobalA.data());
339
340 // Point pt;
341 // for (felInt i = 0; i < localSize; ++i) {
342 // m_linearProblem[1]->findCoordinateWithIdDof(idGlobalA[i],pt);
343 // value[i] = 17827 - FelisceParam::instance().density * 980.665 * (pt.y()+ 2.638);
344 // }
345 // m_linearProblem[1]->solution().setValues(localSize,idGlobalP.data(),value.data(),INSERT_VALUES);
346 // m_linearProblem[1]->solution().assembly();
347 // m_linearProblem[1]->gatherSolution();
348 // }
349
350 // Write solution for postprocessing (if required)
351 100 writeSolution();
352
353 // Advance time step.
354 100 updateTime();
355 // Print time information
356 100 printNewTimeIterationBanner();
357
358
2/2
✓ Branch 1 taken 60 times.
✓ Branch 2 taken 40 times.
 100 if( FelisceParam::instance().orderPressureExtrapolation > 0 )
359 60 m_prevPress.copyFrom(m_linearProblem[1]->solution());
360
361
2/2
✓ Branch 0 taken 200 times.
✓ Branch 1 taken 100 times.
 300 for(std::size_t iLinPb=0; iLinPb<2; iLinPb++) {
362
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
 200 if (FelisceParam::verbose())
363 ✗ PetscPrintf(MpiInfo::petscComm(), "NSFracStep: linear problem #%ld\n",static_cast<unsigned long>(iLinPb));
364
365 // first linear problem
366
2/2
✓ Branch 0 taken 100 times.
✓ Branch 1 taken 100 times.
 200 if ( iLinPb==0 ) {
367
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 92 times.
 100 if ( m_fstransient->iteration == 1 ) {
368 // Initialize of bdf solver.
369 8 m_bdfFS.initialize(m_linearProblem[0]->solution());
370 // Compute the extrapolate.
371 8 m_bdfFS.extrapolate(m_solExtrapolate);
372 // Initialize it.
373 8 m_linearProblem[0]->initExtrapol(m_solExtrapolate);
374
375
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 if(FelisceParam::instance().useALEformulation)
376 4 m_tau = 1./(m_fstransient->timeStep);
377
378
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
 8 if ( FelisceParam::instance().readPreStressedDisplFromFile )
379 ✗ computePreStressedDispMeshVector();
380
381
382 } else {
383 // Update bdf.
384 92 m_bdfFS.update(m_linearProblem[0]->solution());
385 // Compute the extrapolate.
386 92 m_bdfFS.extrapolate(m_solExtrapolate);
387 // Initialize it.
388 92 m_linearProblem[0]->updateExtrapol(m_solExtrapolate);
389 }
390 // convective velocity for supg stabilization in projection step
391 100 m_linearProblem[0]->gatherVector(m_solExtrapolate,m_beta);
392
393 }
394
395 // second problem (projection)
396 // lumpedModelBC bc computation
397
4/4
✓ Branch 0 taken 80 times.
✓ Branch 1 taken 120 times.
✓ Branch 2 taken 40 times.
✓ Branch 3 taken 40 times.
 200 if (m_lumpedModelBC && iLinPb==1) {
398
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 36 times.
 40 if ( m_fstransient->iteration == 1 ) {
399 // m_labelFluxListLumpedModelBC computation with LumpedModelBC labels
400 4 m_labelFluxListLumpedModelBC.clear();
401
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 4 times.
 12 for(std::size_t i=0; i<m_lumpedModelBC->size(); i++) {
402 8 m_labelFluxListLumpedModelBC.push_back(FelisceParam::instance().lumpedModelBCLabel[i]);
403 }
404 }
405 // update the fluxes
406 // justine's warning (10-09-12) : with a velocity with divergence not necesary free... ?
407 40 m_linearProblem[0]->computeMeanQuantity(velocity,m_labelFluxListLumpedModelBC,m_fluxLumpedModelBC);
408
2/2
✓ Branch 1 taken 80 times.
✓ Branch 2 taken 40 times.
 120 for(std::size_t i=0; i<m_lumpedModelBC->size(); i++) {
409 80 m_lumpedModelBC->Q(i) = m_fluxLumpedModelBC[i] ;
410 }
411 // iterate the lumpedModelBC parameters
412 40 m_lumpedModelBC->iterate();
413 }
414
415
2/6
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 200 times.
 200 if (FelisceParam::verbose() and ( iLinPb == 0 ))
416 ✗ PetscPrintf(MpiInfo::petscComm(),"Linear problems NS: Advection-Diffusion-step \n");
417
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
 200 else if (FelisceParam::verbose())
418 ✗ PetscPrintf(MpiInfo::petscComm(),"Linear problems NS: Projection-step \n");
419
420 // Compute the RHS linked to the time scheme (bdf).
421 200 m_bdfFS.computeRHSTime(m_fstransient->timeStep);
422 // Fill _seqBdfRHS.
423 200 m_linearProblem[iLinPb]->gatherVectorBeforeAssembleMatrixRHS();
424
425
6/6
✓ Branch 0 taken 100 times.
✓ Branch 1 taken 100 times.
✓ Branch 3 taken 60 times.
✓ Branch 4 taken 40 times.
✓ Branch 5 taken 60 times.
✓ Branch 6 taken 140 times.
 200 if( iLinPb==0 && FelisceParam::instance().useALEformulation > 0) {
426 60 computeALEterms();
427 } else {
428 //Embedded interface in the projection step
429
4/6
✓ Branch 0 taken 100 times.
✓ Branch 1 taken 40 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 140 times.
 140 if( iLinPb==1 && FelisceParam::instance().EmbeddedInterface){
430 ✗ if ( FelisceParam::instance().fsiCoupling )
431 ✗ m_linearProblem[iLinPb]->meshLocal()->moveMesh(m_dispMeshVector,0.0);
432
433 // Enbedded interface crossed terms are calcualted in the reference configuration
434 ✗ static_cast<LinearProblemNSFracStepProj*>(m_linearProblem[1])->userAssembleMatrixRHSEmbeddedInterface();
435
436 ✗ if ( FelisceParam::instance().fsiCoupling )
437 ✗ m_linearProblem[iLinPb]->meshLocal()->moveMesh(m_dispMeshVector,1.0);
438 }
439
440 //Assembly loop on elements.
441 140 m_linearProblem[iLinPb]->assembleMatrixRHS(MpiInfo::rankProc());
442
443
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 140 times.
 140 if ( FelisceParam::instance().fsiCoupling ) {
444 ✗ if ( iLinPb != 1 ) {
445 ✗ std::cout << "BUG !!!! "<<std::endl;
446 ✗ exit(1);
447 }
448 ✗ m_rhsPressWihoutBC.copyFrom(m_linearProblem[iLinPb]->vector());
449 }
450 }
451
452 //Specific operations before solving the system.
453 200 postAssemblingMatrixRHS(iLinPb);
454
455
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
 200 if ( FelisceParam::instance().fsiCoupling )
456 ✗ m_linearProblem[iLinPb]->meshLocal()->moveMesh(m_dispMeshVector,0.0);
457
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
 200 if ( FelisceParam::instance().readPreStressedDisplFromFile )
458 ✗ m_linearProblem[iLinPb]->meshLocal()->moveMesh(m_preStressedDispMeshVector,1.0);
459
460
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
 200 if (FelisceParam::verbose())
461 ✗ PetscPrintf(MpiInfo::petscComm(),"Applying BC \n");
462
463 //Apply essential transient boundary conditions.
464 200 m_linearProblem[iLinPb]->finalizeEssBCTransient();
465
466 // Apply all the boundary conditions.
467 200 m_linearProblem[iLinPb]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc());
468
469
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
 200 if ( FelisceParam::instance().readPreStressedDisplFromFile )
470 ✗ m_linearProblem[iLinPb]->meshLocal()->moveMesh(m_preStressedDispMeshVector,-1.0);
471
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
 200 if ( FelisceParam::instance().fsiCoupling )
472 ✗ m_linearProblem[iLinPb]->meshLocal()->moveMesh(m_dispMeshVector,1.0);
473
474 //Solve linear system.
475 200 m_linearProblem[iLinPb]->solve(MpiInfo::rankProc(), MpiInfo::numProc());
476
477 //Compute pressure if incremental scheme is used (note that we solve for the increment \delta p )
478
6/6
✓ Branch 0 taken 100 times.
✓ Branch 1 taken 100 times.
✓ Branch 3 taken 60 times.
✓ Branch 4 taken 40 times.
✓ Branch 5 taken 60 times.
✓ Branch 6 taken 140 times.
 200 if ( ( iLinPb == 1 ) && (FelisceParam::instance().orderPressureExtrapolation > 0) ) { // && (m_fstransient->iteration > 1)
479 // Correction of the pressure: p^n = p^{n-1} + \delta p^n
480 60 m_linearProblem[1]->solution().axpy(1, m_prevPress);
481 }
482
483 //Gather solution
484 200 m_linearProblem[iLinPb]->gatherSolution();
485
486 }
487
488
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 100 times.
 100 if ( FelisceParam::instance().fsiCoupling )
489 ✗ computeFluidLoad();
490 100 }
491
492 ✗ int NSFracStepModel::getNstate() const {
493 ✗ return m_linearProblem[0]->numDof();
494 }
495
496 ✗ void NSFracStepModel::getState(double* & state) {
497 ✗ m_linearProblem[0]->getSolution(state, MpiInfo::numProc(), MpiInfo::rankProc());
498 }
499
500 ✗ felInt NSFracStepModel::numDof_swig() {
501 ✗ return m_linearProblem[0]->numDof();
502 }
503
504 ✗ void NSFracStepModel::getState_swig(double* data, felInt numDof) {
505 double * state;
506 ✗ m_linearProblem[0]->getSolution(state, MpiInfo::numProc(), MpiInfo::rankProc());
507 ✗ for (felInt i=0 ; i<numDof ; i++) {
508 ✗ data[i] = state[i];
509 }
510 }
511
512 ✗ void NSFracStepModel::setState(double* & state) {
513 ✗ m_linearProblem[0]->setSolution(state, MpiInfo::numProc(), MpiInfo::rankProc());
514 }
515
516
517 // *********************** ALE ***********************
518 60 void NSFracStepModel::computeALEterms() {
519 // Assembly the RHS (flagMatrixRHS = FlagMatrixRHS::only_rhs) in the previous mesh configuration for the AdvDiff problem:
520
521
2/2
✓ Branch 1 taken 56 times.
✓ Branch 2 taken 4 times.
 60 if ( m_fstransient->iteration>1 ) {
522 56 m_linearProblem[0]->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::only_rhs);
523
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
 56 if ( FelisceParam::instance().fsiCoupling )
524 ✗ m_rhsWihoutBC.copyFrom(m_linearProblem[0]->vector());
525 }
526
527 // Solve the extension problem and move the mesh of _linearProblem[0] and _linearProblem[1]
528 // according to the mesh displacement obtained
529
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if( FelisceParam::instance().updateMeshByVelocity ) {
530 // the HarmonicExtention unknown is velocity
531 // dispalcement is computed analitycally
532 ✗ UpdateMeshByVelocity();
533 } else {
534 // the HarmonicExtention unknown is displacement
535 // velocity is computed analitycally
536 60 UpdateMeshByDisplacement();
537 }
538
539
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 56 times.
 60 if ( m_fstransient->iteration == 1 ) {
540 4 m_linearProblem[0]->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::only_rhs);
541
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( FelisceParam::instance().fsiCoupling )
542 ✗ m_rhsWihoutBC.copyFrom(m_linearProblem[0]->vector());
543 }
544
545 // m_beta = u^{n-1}-w^n
546 60 m_beta.axpy(1.,m_velMesh);
547
548 // Assembly the Matrix (flagMatrixRHS = FlagMatrixRHS::only_matrix) in the current mesh configuration
549 // for the AdvDiff problem
550 60 m_linearProblem[0]->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::only_matrix);
551
552
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if ( FelisceParam::instance().fsiCoupling ) {
553
554 // Saturate the matrix to impose cero Drichlet conditions on \pd \Sigma
555 ✗ if ( m_hasLumpedMassRN && ( m_fstransient->iteration == 1 ) && (FelisceParam::instance().fsiRNscheme != -1) ) {
556 ✗ PetscScalar dti = 1./m_fstransient->timeStep;
557 ✗ m_lumpedMassRN.scale(dti);
558 }
559 ✗ if ( FelisceParam::instance().fsiRNscheme != -1)
560 ✗ m_linearProblem[0]->matrix().diagonalSet(m_lumpedMassRN,ADD_VALUES);
561
562 ✗ if (m_matrixWithoutBCcreated == 0) {
563 ✗ m_matrixWithoutBC.duplicateFrom(m_linearProblem[0]->matrix(),MAT_COPY_VALUES);
564 ✗ m_matrixWithoutBCcreated = true;
565 }
566 ✗ m_matrixWithoutBC.copyFrom(m_linearProblem[0]->matrix(), SAME_NONZERO_PATTERN);
567
568 }
569 60 }
570
571
572 ✗ void NSFracStepModel::UpdateMeshByVelocity() {
573
574 ✗ if(m_fstransient->iteration>1) {
575 ✗ m_solDispl.copyFrom(m_linearProblem[2]->solution());
576 }
577
578 ✗ felInt displIterMax = FelisceParam::instance().timeMaxCaseFile; // to do: it should be taken automatically form displacement file.case
579
580 ✗ if(m_fstransient->iteration == m_cycl * displIterMax) { // to read the file case in cycle
581 ✗ m_displTimeStep = 0.;
582 ✗ m_cycl++;
583 }
584 ✗ PetscPrintf(MpiInfo::petscComm(),"Cycle %d, Read data at time t = %f\n", m_cycl, m_displTimeStep);
585 ✗ static_cast<LinearProblemHarmonicExtension*>(m_linearProblem[2])->readDataDisplacement(m_ios, m_displTimeStep); // defined in linearProblemHarmonicExtention
586
587 //Solve the extension problem
588 ✗ solveHarmExt();
589
590 // Store velocity and apply it to the NS problem. seqSolution is the velocity
591 ✗ m_linearProblem[2]->gatherSolution();
592 ✗ m_velMesh.copyFrom(m_linearProblem[2]->sequentialSolution());
593 ✗ m_velMesh.scale(-1.);
594 //Compute displacement: d^n = w^n*dt + d^{n-1}
595 //VecAXPBY(Vec y,PetscScalar a,PetscScalar b,Vec x); -> y=a∗x+b∗y
596 // solution() == w^n
597 //_solDispl == d^{n-1}
598 ✗ m_linearProblem[2]->solution().axpby(1, FelisceParam::instance().timeStep, m_solDispl);
599 ✗ m_linearProblem[2]->gatherSolution();
600 // seqSolution is now the displacement
601
602 ✗ updateMesh();
603 }
604
605 60 void NSFracStepModel::UpdateMeshByDisplacement() {
606
607 // fsi challenge
608
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if ( FelisceParam::instance().useElasticExtension )
609 ✗ m_dispMesho_pp.copyFrom( m_linearProblem[2]->solution() );
610
611 // m_dispMesh = 0 at m_fstransient->iteration = 1
612 60 m_velMesh.copyFrom(m_dispMesh);
613 //Solve the extension problem
614 60 solveHarmExt();
615 //Gather solution
616 60 m_linearProblem[2]->gatherSolution();
617 //Update mesh
618 60 updateMesh();
619
620 // Compute velocity and apply it to the NS problem. seqSolution is the displacement
621 // VecAXPBY(Vec y,PetscScalar a,PetscScalar b,Vec x); -> y=a∗x+b∗y
622 // now m_velMesh == d^{n-1}
623 // and m_dispMesh == d^{n}
624
2/2
✓ Branch 1 taken 56 times.
✓ Branch 2 taken 4 times.
 60 if ( m_fstransient->iteration > 1 )
625 56 m_velMesh.axpby(-m_tau,m_tau,m_dispMesh);
626 // Result: m_velMesh == -w^{n} = \dfrac{1}{\dt} (-d^{n} + d^{n-1})
627 //double norm;
628 //VecNorm(m_velMesh.toPetsc(), NORM_INFINITY ,&norm);
629 //std::cout << std::endl << " Norm vel mesh = " << norm << std::endl <<std::endl;
630
631 60 }
632
633
634 60 void NSFracStepModel::solveHarmExt() {
635 //Solve the extension problem
636 //Assemble matrix and RHS
637 60 m_linearProblem[2]->assembleMatrixRHS(MpiInfo::rankProc());
638 //Specific operations before solve the system
639 60 postAssemblingMatrixRHS(2);
640 //Apply essential transient boundary conditions
641 60 m_linearProblem[2]->finalizeEssBCTransient();
642 60 m_linearProblem[2]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc());
643 //Solve linear system
644 60 m_linearProblem[2]->solve(MpiInfo::rankProc(), MpiInfo::numProc());
645
646 // fsi challenge (update displacement with increment)
647
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if ( FelisceParam::instance().useElasticExtension )
648 ✗ m_linearProblem[2]->solution().axpy(1.,m_dispMesho_pp);
649
650 60 }
651
652 60 void NSFracStepModel::updateMesh() {
653
654
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if (FelisceParam::verbose())
655 ✗ PetscPrintf(MpiInfo::petscComm(),"Update mesh\n");
656
657 60 m_dispMesh.copyFrom(m_linearProblem[2]->sequentialSolution());
658 // _dispMesh is the solution obtained with solveHarmExt()
659 // Extract the values of _dispMesh and put them in the std::vector _dispMeshVector
660 // according to the application ordering of the nodes
661 60 m_dispMesh.getValues(m_numDofExtensionProblem, m_petscToGlobal.data(), m_tmpDisp.data());
662
663 60 felInt numNodes = m_linearProblem[2]->supportDofUnknown(0).listNode().size();
664 felInt idof;
665 60 felInt numCompDisp = m_linearProblem[2]->dimension();
666 (void)numCompDisp;
667 60 felInt numCoorMesh = m_linearProblem[2]->meshLocal()->numCoor();
668
669
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if(FelisceParam::instance().FusionDof) {
670 ✗ for (felInt ip = 0; ip < numNodes; ip++) {
671 ✗ idof = m_linearProblem[1]->supportDofUnknown(0).listPerm()[ip];
672 ✗ if ( numCompDisp < numCoorMesh ) {
673 ✗ m_dispMeshVector[ip*3] = m_tmpDisp[idof*2];
674 ✗ m_dispMeshVector[ip*3+1] = m_tmpDisp[idof*2+1];
675 ✗ m_dispMeshVector[ip*3+2] = 0.;
676 }
677 ✗ else if ( numCompDisp == numCoorMesh )
678 ✗ for (felInt ic=0; ic < numCoorMesh; ++ic)
679 ✗ m_dispMeshVector[ip*numCoorMesh+ic] = m_tmpDisp[idof*numCoorMesh+ic];
680 }
681 }
682 else
683 60 m_dispMeshVector = m_tmpDisp;
684
685 //Move the mesh in the AdvDiff and Proj problems
686
1/2
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
 60 m_linearProblem[0]->meshLocal()->moveMesh(m_dispMeshVector,1.0);
687
1/2
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
 60 m_linearProblem[1]->meshLocal()->moveMesh(m_dispMeshVector,1.0);
688
689
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
 60 if ( FelisceParam::instance().useElasticExtension )
690 // extension computed on the deformed mesh (takes into account the size of the deformed elements)
691 ✗ m_linearProblem[2]->meshLocal()->moveMesh(m_dispMeshVector,1.0);
692
693 60 }
694
695
696 // *********************** COMPUTE PRE-STRESSED DISP ***********************
697 ✗ void NSFracStepModel::computePreStressedDispMeshVector() {
698
699 ✗ static_cast<LinearProblemHarmonicExtension*>(m_linearProblem[2])->HarmExtSol().getValues(m_numDofExtensionProblem, m_petscToGlobal.data(), m_tmpDisp.data());
700
701 ✗ felInt numNodes = m_linearProblem[2]->supportDofUnknown(0).listNode().size();
702 ✗ if(FelisceParam::instance().FusionDof) {
703 // if there are e.g. RIS model, the number of dofs are re-numerated.
704 ✗ m_preStressedDispMeshVector.resize(numNodes*3, 0.);
705 ✗ for (felInt ip = 0; ip < numNodes; ip++) {
706 ✗ m_preStressedDispMeshVector[ip*3] = m_tmpDisp[m_linearProblem[2]->supportDofUnknown(0).listPerm()[ip]*3];
707 ✗ m_preStressedDispMeshVector[ip*3+1] = m_tmpDisp[m_linearProblem[2]->supportDofUnknown(0).listPerm()[ip]*3+1];
708 ✗ m_preStressedDispMeshVector[ip*3+2] = m_tmpDisp[m_linearProblem[2]->supportDofUnknown(0).listPerm()[ip]*3+2];
709 }
710 } else
711 ✗ m_preStressedDispMeshVector = m_tmpDisp;
712 }
713
714 // *********************** FSI ***********************
715 ✗ void NSFracStepModel::computeFluidLoad() {
716
717 // Viscous residual calculation. We calculate -ViscousResidual
718 ✗ m_matrixWithoutBC.scale(-1.);
719 ✗ multAdd(m_matrixWithoutBC,m_linearProblem[0]->solution(),m_rhsWihoutBC,m_viscousResidual);
720 ✗ m_linearProblem[0]->gatherVector(m_viscousResidual, m_seqViscousResidual);
721
722 // Pressure load calculation in FWI fashion
723 ✗ static_cast<LinearProblemNSFracStepAdvDiff*>(m_linearProblem[0])->pressureFWI() = 1;
724 ✗ m_linearProblem[0]->vector().set(0.);
725 ✗ m_linearProblem[0]->assembleMatrixRHSNaturalBoundaryCondition(MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, 0);
726 //_linearProblem[0]->gatherVector(_linearProblem[0]->RHS(), _seqPressLoad);
727 ✗ static_cast<LinearProblemNSFracStepAdvDiff*>(m_linearProblem[0])->pressureFWI() = 0;
728
729 // PressureFWI - ViscousResidual
730 ✗ waxpy(m_residual,1., m_linearProblem[0]->vector(),m_viscousResidual);
731 ✗ m_linearProblem[0]->gatherVector(m_residual, m_seqResidual);
732
733 }
734
735 ✗ PetscVector& NSFracStepModel::seqResidual() {
736 ✗ return m_seqResidual;
737 }
738 ✗ PetscVector& NSFracStepModel::residual() {
739 ✗ return m_residual;
740 }
741 ✗ PetscVector& NSFracStepModel::seqViscousResidual() {
742 ✗ return m_seqViscousResidual;
743 }
744
745 ✗ int& NSFracStepModel::hasLumpedMassRN() {
746 ✗ return m_hasLumpedMassRN;
747 }
748 ✗ PetscVector& NSFracStepModel::lumpedMassRN() {
749 ✗ return m_lumpedMassRN;
750 }
751
752 // *********************** Incremental scheme ***********************
753 ✗ void NSFracStepModel::computePressure(){
754 // Correction of the pressure: p^n = p^{n-1} + \delta p^n
755 ✗ m_linearProblem[1]->solution().axpy(1, m_prevPress);
756 ✗ m_prevPress.copyFrom(m_linearProblem[1]->solution());
757 }
758
759 }
760