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: Julien Castelneau, Dominique Chapelle, Miguel Fernandez |
13 |
|
|
// |
14 |
|
|
|
15 |
|
|
// System includes |
16 |
|
|
|
17 |
|
|
// External includes |
18 |
|
|
|
19 |
|
|
// Project includes |
20 |
|
|
#include "NSContinuationModel.hpp" |
21 |
|
|
|
22 |
|
|
|
23 |
|
|
namespace felisce { |
24 |
|
✗ |
NSContinuationModel::NSContinuationModel():Model() { |
25 |
|
✗ |
m_name = "NSContinuation"; |
26 |
|
|
} |
27 |
|
|
|
28 |
|
|
// Called in Model::initializeLinearProblem, just after the mesh partitioning. |
29 |
|
✗ |
void NSContinuationModel::initializeDerivedModel() { |
30 |
|
|
// get the linearProblem for continuation; m_linearProblem[0] = continuation and m_linearProblem[1] = harmonic extension |
31 |
|
✗ |
m_lpb = static_cast<LinearProblemNSContinuation*>(m_linearProblem[0]); |
32 |
|
|
|
33 |
|
✗ |
m_tau = 1./(m_fstransient->timeStep); |
34 |
|
|
|
35 |
|
|
// Initializing the PetscVectors |
36 |
|
✗ |
m_lpb->initPetscVectors(); |
37 |
|
|
|
38 |
|
|
// Initialize the displacement vector |
39 |
|
✗ |
if( m_lpb->dimension()==2) |
40 |
|
|
{ |
41 |
|
✗ |
m_disp.duplicateFrom(m_lpb->vector()); |
42 |
|
✗ |
m_disp.zeroEntries(); |
43 |
|
|
} |
44 |
|
|
|
45 |
|
|
|
46 |
|
✗ |
if( FelisceParam::instance(instanceIndex()).useALEformulation ) |
47 |
|
|
{ |
48 |
|
|
// ALE |
49 |
|
✗ |
m_numDofExtensionProblem = m_linearProblem[1]->numDof(); |
50 |
|
✗ |
m_dispMeshVector.resize( m_linearProblem[1]->supportDofUnknown(0).listNode().size() * m_linearProblem[1]->meshLocal()->numCoor() ); |
51 |
|
✗ |
m_petscToGlobal.resize(m_numDofExtensionProblem); |
52 |
|
✗ |
for (int ipg=0; ipg< m_numDofExtensionProblem; ipg++) |
53 |
|
✗ |
m_petscToGlobal[ipg] = ipg; |
54 |
|
|
|
55 |
|
✗ |
AOApplicationToPetsc(m_linearProblem[1]->ao(),m_numDofExtensionProblem,m_petscToGlobal.data()); |
56 |
|
|
} |
57 |
|
|
|
58 |
|
|
|
59 |
|
|
} |
60 |
|
|
|
61 |
|
✗ |
void NSContinuationModel::setExternalVec() { |
62 |
|
|
|
63 |
|
|
// ALE |
64 |
|
✗ |
if( FelisceParam::instance(instanceIndex()).useALEformulation ) |
65 |
|
|
{ |
66 |
|
|
// Initialization of the std::vector containing the displacement of the mesh d^{n}, sequential std::vector |
67 |
|
✗ |
m_dispMesh.duplicateFrom(m_linearProblem[1]->sequentialSolution()); |
68 |
|
✗ |
m_dispMesh.zeroEntries(); |
69 |
|
|
|
70 |
|
|
// Initialization of the std::vector containing the velocity of the mesh w, sequential std::vector |
71 |
|
✗ |
m_velMesh.duplicateFrom(m_linearProblem[1]->sequentialSolution()); |
72 |
|
✗ |
m_velMesh.zeroEntries(); |
73 |
|
|
|
74 |
|
|
// Passing the ao and dof |
75 |
|
✗ |
m_lpb->pushBackExternalAO(m_linearProblem[1]->ao()); |
76 |
|
✗ |
m_lpb->pushBackExternalDof(m_linearProblem[1]->dof()); |
77 |
|
|
// Passing the velocity of the mesh to the NS problem |
78 |
|
✗ |
m_lpb->pushBackExternalVec(m_velMesh); //externalVec(0) in m_lpb |
79 |
|
|
// Resizing the auxiliary vectors needed in updateMesh() |
80 |
|
✗ |
m_tmpDisp.resize(m_numDofExtensionProblem); |
81 |
|
|
|
82 |
|
✗ |
m_displTimeStep = 0.; |
83 |
|
|
|
84 |
|
|
} |
85 |
|
|
|
86 |
|
✗ |
if(m_lpb->dimension()==2) |
87 |
|
|
{ |
88 |
|
✗ |
m_dispSeq.createSeq(PETSC_COMM_SELF, m_lpb->numDof()); |
89 |
|
✗ |
m_dispSeq.set(0.); |
90 |
|
✗ |
m_lpb->pushBackExternalVec(m_dispSeq); //externalVec(1) if in 2D |
91 |
|
|
} |
92 |
|
|
|
93 |
|
|
} |
94 |
|
|
|
95 |
|
|
|
96 |
|
✗ |
void NSContinuationModel::preAssemblingMatrixRHS(std::size_t iProblem) { |
97 |
|
✗ |
if(iProblem == 0) { |
98 |
|
✗ |
m_linearProblem[0]->solution().zeroEntries(); |
99 |
|
|
} |
100 |
|
|
|
101 |
|
|
//========= |
102 |
|
|
// ALE |
103 |
|
|
//========= |
104 |
|
✗ |
if( FelisceParam::instance(instanceIndex()).useALEformulation ) { |
105 |
|
✗ |
if ( ( iProblem == 1 ) && ( !FelisceParam::instance(instanceIndex()).useElasticExtension ) ) { |
106 |
|
|
// First assembly loop in iteration 0 to build static matrix of the extension problem |
107 |
|
✗ |
m_linearProblem[1]->assembleMatrixRHS(MpiInfo::rankProc()); |
108 |
|
|
// Save the static part of the extension problem in the matrix _A |
109 |
|
✗ |
m_linearProblem[1]->copyMatrixRHS(); |
110 |
|
|
} |
111 |
|
|
} |
112 |
|
|
} |
113 |
|
|
|
114 |
|
|
|
115 |
|
|
//! The main purpose of this function is to sum the constant and the non constant matrices |
116 |
|
✗ |
void NSContinuationModel::postAssemblingMatrixRHS(std::size_t iProblem) { |
117 |
|
|
// matrix(0) = matrix(0) + matrix(1) (add static matrix to the dynamic matrix to build complete matrix of the system ) |
118 |
|
✗ |
m_linearProblem[iProblem]->addMatrixRHS(); |
119 |
|
|
} |
120 |
|
|
|
121 |
|
|
|
122 |
|
✗ |
void NSContinuationModel::prepareForward() { |
123 |
|
|
|
124 |
|
✗ |
writeSolution(); |
125 |
|
|
|
126 |
|
✗ |
updateTime(); //Here the matrix(0), the non constant one, is cleared and the rhs also. |
127 |
|
|
|
128 |
|
✗ |
printNewTimeIterationBanner(); |
129 |
|
|
|
130 |
|
✗ |
m_lpb->gatherVectorBeforeAssembleMatrixRHS(); |
131 |
|
|
|
132 |
|
✗ |
m_lpb->readVelocityData(*io(), m_fstransient->time); |
133 |
|
|
|
134 |
|
✗ |
computeALEterms(); // if no ALE, it does not compute the displacement |
135 |
|
|
} |
136 |
|
|
|
137 |
|
|
|
138 |
|
✗ |
void NSContinuationModel::solveContinuation() { |
139 |
|
✗ |
const auto& r_felisce_param_instance = FelisceParam::instance(instanceIndex()); |
140 |
|
|
|
141 |
|
✗ |
m_lpb->solution().zeroEntries(); |
142 |
|
|
|
143 |
|
|
//Steps corresponding to solveNS() |
144 |
|
✗ |
if( FelisceParam::instance(instanceIndex()).useALEformulation ) |
145 |
|
✗ |
m_lpb->meshLocal()->moveMesh(m_dispMeshVector,0.0); |
146 |
|
|
|
147 |
|
|
//Apply boundary conditions for the continuation problem in the reference mesh configuration |
148 |
|
✗ |
m_lpb->finalizeEssBCTransient(); |
149 |
|
✗ |
m_lpb->applyBC(r_felisce_param_instance.essentialBoundaryConditionsMethod, MpiInfo::rankProc()); |
150 |
|
|
|
151 |
|
✗ |
if( FelisceParam::instance(instanceIndex()).useALEformulation ) |
152 |
|
✗ |
m_lpb->meshLocal()->moveMesh(m_dispMeshVector,1.0); |
153 |
|
|
|
154 |
|
|
// Solve the linear system |
155 |
|
✗ |
m_lpb->solve(MpiInfo::rankProc(), MpiInfo::numProc()); |
156 |
|
|
|
157 |
|
|
/*if(m_lpb->dimension() == 2) |
158 |
|
|
{ |
159 |
|
|
m_disp.axpy(m_fstransient->timeStep, m_lpb->solution()); |
160 |
|
|
m_lpb->gatherVector(m_disp, m_dispSeq); |
161 |
|
|
} */ |
162 |
|
|
|
163 |
|
|
// Gather the velocity |
164 |
|
✗ |
m_lpb->gatherSolution(); |
165 |
|
|
} |
166 |
|
|
|
167 |
|
✗ |
void NSContinuationModel::forward() { |
168 |
|
✗ |
prepareForward(); |
169 |
|
✗ |
solveContinuation(); |
170 |
|
|
} |
171 |
|
|
|
172 |
|
|
// *********************** ALE ************************************************* |
173 |
|
✗ |
void NSContinuationModel::computeALEterms() { |
174 |
|
|
// Assemble the RHS in the previous mesh configuration |
175 |
|
✗ |
m_lpb->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::only_rhs); |
176 |
|
|
|
177 |
|
✗ |
if( FelisceParam::instance(instanceIndex()).useALEformulation ) |
178 |
|
✗ |
UpdateMeshByDisplacement(); |
179 |
|
|
|
180 |
|
✗ |
m_lpb->assembleCIPStabilization(); |
181 |
|
|
|
182 |
|
|
// Assemble the matrix in the current mesh configuration |
183 |
|
✗ |
m_lpb->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::matrix_and_rhs); |
184 |
|
|
} |
185 |
|
|
|
186 |
|
✗ |
void NSContinuationModel::UpdateMeshByDisplacement() { |
187 |
|
|
// Get the previous displacement |
188 |
|
✗ |
m_velMesh.copyFrom(m_dispMesh); |
189 |
|
|
|
190 |
|
✗ |
if( FelisceParam::instance(instanceIndex()).useALEformulation ) { |
191 |
|
|
// Solve the harmonic extension problem |
192 |
|
✗ |
solveHarmExt(); |
193 |
|
✗ |
m_linearProblem[1]->gatherSolution(); |
194 |
|
|
} |
195 |
|
✗ |
updateMesh(); // current displacement m_dispMesh |
196 |
|
|
|
197 |
|
|
// Update the mesh velocity; VecAXPBY(Vec y,PetscScalar a,PetscScalar b,Vec x); -> y=a∗x+b∗y |
198 |
|
✗ |
m_velMesh.axpby(-m_tau,m_tau,m_dispMesh); |
199 |
|
|
// Result: m_velMesh = -w^{n} = \dfrac{1}{\dt} (-d^{n} + d^{n-1}) |
200 |
|
|
} |
201 |
|
|
|
202 |
|
✗ |
void NSContinuationModel::solveHarmExt() { |
203 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Solve harmonic extension\n"); |
204 |
|
|
//Assemble matrix and RHS |
205 |
|
✗ |
m_linearProblem[1]->assembleMatrixRHS(MpiInfo::rankProc()); |
206 |
|
|
|
207 |
|
|
//Specific operations before solving the system |
208 |
|
✗ |
postAssemblingMatrixRHS(1); |
209 |
|
|
|
210 |
|
|
//Apply essential transient boundary conditions |
211 |
|
✗ |
m_linearProblem[1]->finalizeEssBCTransient(); |
212 |
|
✗ |
m_linearProblem[1]->applyBC(FelisceParam::instance(instanceIndex()).essentialBoundaryConditionsMethod, MpiInfo::rankProc()); |
213 |
|
|
|
214 |
|
|
//Solve the linear system |
215 |
|
✗ |
m_linearProblem[1]->solve(MpiInfo::rankProc(), MpiInfo::numProc()); |
216 |
|
|
|
217 |
|
|
} |
218 |
|
|
|
219 |
|
✗ |
void NSContinuationModel::updateMesh() { |
220 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Update mesh\n"); |
221 |
|
|
|
222 |
|
|
// Get the displacement from the harmonic extension |
223 |
|
✗ |
m_dispMesh.copyFrom(m_linearProblem[1]->sequentialSolution()); |
224 |
|
|
|
225 |
|
|
// Extract the values of m_dispMesh and put them in the std::vector m_dispMeshVector according to the application ordering of the nodes |
226 |
|
✗ |
m_dispMesh.getValues(m_numDofExtensionProblem, m_petscToGlobal.data(), m_tmpDisp.data()); |
227 |
|
|
|
228 |
|
|
|
229 |
|
✗ |
m_dispMeshVector = m_tmpDisp; |
230 |
|
|
|
231 |
|
|
//Move the mesh in the fluid problem |
232 |
|
✗ |
m_lpb->meshLocal()->moveMesh(m_dispMeshVector,1.0); |
233 |
|
|
} |
234 |
|
|
|
235 |
|
|
// Compute the L2-error and append to file [Temporary implementation] |
236 |
|
|
|
237 |
|
|
|
238 |
|
|
} |
239 |
|
|
|