| Directory: | ./ |
|---|---|
| File: | Model/NSModel.cpp |
| Date: | 2024-04-14 07:32:34 |
| Exec | Total | Coverage | |
|---|---|---|---|
| Lines: | 341 | 603 | 56.6% |
| Branches: | 225 | 544 | 41.4% |
| 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 "Core/felisceTools.hpp" | ||
| 21 | #include "Model/NSModel.hpp" | ||
| 22 | #include "Solver/CVGraphInterface.hpp" | ||
| 23 | #include "Solver/linearProblemNS.hpp" | ||
| 24 | #include "Solver/linearProblemHarmonicExtension.hpp" | ||
| 25 | #ifdef FELISCE_WITH_CVGRAPH | ||
| 26 | #include "Model/cvgMainSlave.hpp" | ||
| 27 | #endif | ||
| 28 | |||
| 29 | namespace felisce { | ||
| 30 | |||
| 31 |
17/34✓ Branch 2 taken 42 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 42 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 42 times.
✗ Branch 9 not taken.
✓ Branch 13 taken 42 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 42 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 42 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 42 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 42 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 42 times.
✗ Branch 29 not taken.
✓ Branch 33 taken 42 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 42 times.
✗ Branch 37 not taken.
✓ Branch 39 taken 42 times.
✗ Branch 40 not taken.
✓ Branch 46 taken 42 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 42 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 42 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 42 times.
✗ Branch 56 not taken.
✓ Branch 58 taken 42 times.
✗ Branch 59 not taken.
|
42 | NSModel::NSModel():Model() |
| 32 | { | ||
| 33 | 42 | m_lumpedModelBC = nullptr; | |
| 34 | 42 | m_risModel = nullptr; | |
| 35 |
1/2✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
|
42 | m_name = "NS"; |
| 36 | 42 | m_matrixWithoutBCcreated = false; | |
| 37 | 42 | m_CrankNicolsonNS =nullptr; | |
| 38 | 42 | } | |
| 39 | |||
| 40 | /***********************************************************************************/ | ||
| 41 | /***********************************************************************************/ | ||
| 42 | |||
| 43 | 116 | NSModel::~NSModel() | |
| 44 | { | ||
| 45 | 84 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 46 | |||
| 47 | // Write displacement and normal | ||
| 48 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
|
84 | if ( r_instance.useFDlagrangeMult ) { |
| 49 | ✗ | const int ivar = m_pSolverNS->listVariable().getVariableIdList(lagMultiplier); | |
| 50 | ✗ | const int imesh = m_pSolverNS->listVariable()[ivar].idMesh(); | |
| 51 | ✗ | writeCurrentDisplacement(imesh, "displacement"); | |
| 52 | ✗ | writeCurrentNormalVector(); | |
| 53 | |||
| 54 | ✗ | if ( r_instance.useFicItf ) | |
| 55 | ✗ | writeCurrentDisplacement(2, "dispficstruc"); // TODO D.C. | |
| 56 | } | ||
| 57 | |||
| 58 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 38 times.
|
84 | if ( m_lumpedModelBC ) |
| 59 |
1/2✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
|
8 | delete m_lumpedModelBC; |
| 60 | 116 | } | |
| 61 | |||
| 62 | /***********************************************************************************/ | ||
| 63 | /***********************************************************************************/ | ||
| 64 | |||
| 65 | 42 | void NSModel::setInitialCondition() | |
| 66 | { | ||
| 67 | // Retrieve instance | ||
| 68 | 42 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 69 | |||
| 70 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
42 | if ( r_instance.hasInitialCondition && r_instance.useALEformulation ) { |
| 71 | |||
| 72 | // TODO D.C. only in case of ALE we can have an initial condition??? i don't think so... | ||
| 73 | ✗ | int ivar = m_pSolverHE->listVariable().getVariableIdList(displacement); | |
| 74 | ✗ | int iunk = m_pSolverHE->listUnknown().getUnknownIdList(displacement); | |
| 75 | ✗ | int imesh = m_pSolverHE->listVariable()[ivar].idMesh(); | |
| 76 | |||
| 77 | // displacement | ||
| 78 | int numDof; | ||
| 79 | ✗ | if( r_instance.FusionDof ) | |
| 80 | ✗ | numDof = m_pSolverHE->supportDofUnknown(iunk).listNode().size() * m_pSolverHE->meshLocal(imesh)->numCoor(); | |
| 81 | else | ||
| 82 | ✗ | numDof = m_pSolverHE->numDof(); | |
| 83 | |||
| 84 | // displacement | ||
| 85 | ✗ | std::vector<double> vector( numDof, 0.); | |
| 86 | ✗ | m_ios[imesh]->readVariable(2, 0.0, vector.data(), numDof); | |
| 87 | |||
| 88 | ✗ | if( r_instance.FusionDof ) { | |
| 89 | ✗ | int numCompDisp = m_pSolverHE->dimension(); | |
| 90 | ✗ | int numCoorMesh = m_pSolverHE->meshLocal(imesh)->numCoor(); | |
| 91 | ✗ | int numNodes = m_pSolverHE->supportDofUnknown(iunk).listNode().size(); | |
| 92 | int idof; | ||
| 93 | ✗ | std::vector<double> aux(vector); | |
| 94 | ✗ | for (int ip = 0; ip < numNodes; ip++) { | |
| 95 | ✗ | idof = m_pSolverHE->supportDofUnknown(iunk).listPerm()[ip]; | |
| 96 | ✗ | for (felInt ic=0; ic < numCompDisp; ++ic) | |
| 97 | ✗ | vector[idof*numCompDisp+ic] = aux[ip*numCoorMesh+ic]; | |
| 98 | } | ||
| 99 | } | ||
| 100 | |||
| 101 | int low, high; | ||
| 102 | ✗ | VecGetOwnershipRange(m_pSolverHE->solution().toPetsc(),&low,&high); | |
| 103 | |||
| 104 | ✗ | int localSize = high-low; | |
| 105 | ✗ | std::vector<int> ix(localSize); | |
| 106 | ✗ | std::vector<double> vv(localSize); | |
| 107 | ✗ | std::iota(ix.begin(),ix.end(),low); // fills ix with values from low to high | |
| 108 | |||
| 109 | ✗ | AOPetscToApplication(m_pSolverHE->ao(),localSize,ix.data()); | |
| 110 | ✗ | for (int i=0; i<localSize;++i) | |
| 111 | ✗ | vv[i] = vector[ix[i]]; | |
| 112 | ✗ | AOApplicationToPetsc(m_pSolverHE->ao(),localSize,ix.data()); | |
| 113 | |||
| 114 | ✗ | m_pSolverHE->solution().setValues(localSize,ix.data(),vv.data(),INSERT_VALUES); | |
| 115 | ✗ | m_pSolverHE->solution().assembly(); | |
| 116 | ✗ | m_pSolverHE->gatherSolution(); | |
| 117 | ✗ | updateMesh(); | |
| 118 | } | ||
| 119 | 42 | } | |
| 120 | |||
| 121 | /***********************************************************************************/ | ||
| 122 | /***********************************************************************************/ | ||
| 123 | |||
| 124 | 42 | void NSModel::initializeDerivedModel() | |
| 125 | { | ||
| 126 | // Retrieve instance | ||
| 127 | 42 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 128 | |||
| 129 | 42 | m_pSolverNS = static_cast<LinearProblemNS*>(m_linearProblem[0]); | |
| 130 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 24 times.
|
42 | if ( r_instance.useALEformulation ) |
| 131 | 18 | m_pSolverHE = static_cast<LinearProblemHarmonicExtension*>(m_linearProblem[1]); | |
| 132 | |||
| 133 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
|
42 | if (m_initializeCVGraphInterface) { |
| 134 | ✗ | m_cvgraphInterf = felisce::make_shared<CVGraphInterface>(); | |
| 135 | ✗ | m_cvgraphInterf->initialize(&m_pSolverNS->listVariable(), r_instance.CVGraphInterface); | |
| 136 | ✗ | m_cvgraphInterf->identifyIdDof(m_pSolverNS->dof()); | |
| 137 | ✗ | m_pSolverNS->cvgraphInterf() = m_cvgraphInterf; | |
| 138 | } | ||
| 139 | |||
| 140 |
2/2✓ Branch 1 taken 4 times.
✓ Branch 2 taken 38 times.
|
42 | if ( r_instance.lumpedModelBCLabel.size() ) { |
| 141 |
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); |
| 142 | 4 | m_pSolverNS->initializeLumpedModelBC(m_lumpedModelBC); | |
| 143 | } | ||
| 144 | |||
| 145 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
|
42 | if ( r_instance.RISModels ) { |
| 146 | ✗ | m_risModel = new RISModel (m_fstransient, m_linearProblem) ; | |
| 147 | ✗ | m_pSolverNS->initializeRISModel(m_risModel); | |
| 148 | ✗ | m_chek_use_RISforward = false; | |
| 149 | } | ||
| 150 | |||
| 151 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
|
42 | if ( r_instance.CardiacCycle ) { |
| 152 | ✗ | m_cardiacCycle = new CardiacCycle (m_fstransient, m_risModel, m_linearProblem) ; | |
| 153 | ✗ | m_pSolverNS->initializeCardiacCycle(m_cardiacCycle); | |
| 154 | } | ||
| 155 | |||
| 156 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 34 times.
|
42 | if ( r_instance.hasImmersedData ) { |
| 157 |
1/2✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
8 | m_immersed.initialize(m_pSolverNS->mesh()); |
| 158 | } | ||
| 159 | |||
| 160 | // monolithic fsi | ||
| 161 | // m_disp.duplicateFrom(m_pSolverNS->vector()); | ||
| 162 | // m_disp.zeroEntries(); | ||
| 163 | |||
| 164 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 24 times.
|
42 | if ( r_instance.useALEformulation ) { |
| 165 | 18 | const int ivar = m_pSolverHE->listVariable().getVariableIdList(displacement); | |
| 166 | 18 | const int iunk = m_pSolverHE->listUnknown().getUnknownIdList(displacement); | |
| 167 | 18 | const int iMesh = m_pSolverHE->listVariable()[ivar].idMesh(); | |
| 168 | |||
| 169 | 18 | m_numDofExtensionProblem = m_pSolverHE->numDof(); | |
| 170 |
1/2✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
|
18 | m_dispMeshVector.resize( m_pSolverHE->supportDofUnknown(iunk).listNode().size() * m_pSolverHE->meshLocal(iMesh)->numCoor() ); |
| 171 | 18 | m_petscToGlobal.resize(m_numDofExtensionProblem); | |
| 172 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
|
18 | if ( r_instance.readPreStressedDisplFromFile ) |
| 173 | ✗ | m_preStressedDispMeshVector.resize(m_numDofExtensionProblem); | |
| 174 |
2/2✓ Branch 0 taken 252440 times.
✓ Branch 1 taken 18 times.
|
252458 | for (int ipg=0; ipg< m_numDofExtensionProblem; ipg++) |
| 175 | 252440 | m_petscToGlobal[ipg] = ipg; | |
| 176 | |||
| 177 | 18 | AOApplicationToPetsc(m_pSolverHE->ao(),m_numDofExtensionProblem,m_petscToGlobal.data()); | |
| 178 | |||
| 179 |
2/4✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
|
18 | if ( r_instance.readDisplFromFile || r_instance.readPreStressedDisplFromFile ) { |
| 180 | // defined in linearProblemHarmonicExtention | ||
| 181 | ✗ | m_pSolverHE->readMatchFile_ALE(m_pSolverHE->listVariable(), r_instance.MoveMeshMatchFile); | |
| 182 | ✗ | m_pSolverHE->identifyIdDofToMatch_ALE(m_pSolverHE->dof()); | |
| 183 | |||
| 184 | ✗ | if (r_instance.hasInitialCondition || r_instance.readPreStressedDisplFromFile ) { | |
| 185 | ✗ | PetscPrintf(MpiInfo::petscComm(),"\nCycle 1, Read data displacement at time t = 0.00000\n"); | |
| 186 | ✗ | m_pSolverHE->readDataDisplacement(m_ios, 0.0); | |
| 187 | } | ||
| 188 | } | ||
| 189 | } | ||
| 190 | |||
| 191 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
|
42 | if ( r_instance.orderBdfNS > 2 ) |
| 192 | ✗ | FEL_ERROR("BDF not yet implemented for order greater than 2 with Navier-Stokes."); | |
| 193 | |||
| 194 | 42 | m_bdfNS.defineOrder(r_instance.orderBdfNS); | |
| 195 | 42 | m_pSolverNS->initializeTimeScheme(&m_bdfNS); | |
| 196 | |||
| 197 | // Crank-Nicolson time marching | ||
| 198 | // Manage midpoint approximation for the advection term | ||
| 199 | // TODO D.C. does it work??? | ||
| 200 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 42 times.
|
42 | if ( Tools::equal(r_instance.theta,0.5) ) { |
| 201 | ✗ | m_CrankNicolsonNS = new Bdf; | |
| 202 | ✗ | m_CrankNicolsonNS->defineOrder(2); | |
| 203 | // Second order approximation of the midpoint | ||
| 204 | ✗ | m_CrankNicolsonNS->beta(0) = 3./2; | |
| 205 | ✗ | m_CrankNicolsonNS->beta(1) = -1./2; | |
| 206 | } | ||
| 207 | |||
| 208 | // Initialization of residual stuff | ||
| 209 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
42 | if ( r_instance.withCVG && !r_instance.fsiCoupling ) { |
| 210 | ✗ | m_rhsWihoutBC.duplicateFrom(m_pSolverNS->vector()); | |
| 211 | ✗ | m_residual.duplicateFrom(m_pSolverNS->vector()); | |
| 212 | ✗ | m_seqResidual.duplicateFrom(m_pSolverNS->sequentialSolution()); | |
| 213 | } | ||
| 214 | |||
| 215 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 36 times.
|
42 | if ( r_instance.fsiCoupling ) { |
| 216 | // TODO: we can use the functions createAndCopyMatrixRHSWithoutBC and computeResidual in LinearProblem | ||
| 217 | 6 | m_rhsWihoutBC.duplicateFrom(m_pSolverNS->vector()); | |
| 218 | 6 | m_residual.duplicateFrom(m_pSolverNS->vector()); | |
| 219 | 6 | m_lumpedMassRN.duplicateFrom(m_pSolverNS->vector()); | |
| 220 | 6 | m_seqResidual.duplicateFrom(m_pSolverNS->sequentialSolution()); | |
| 221 | } | ||
| 222 | 42 | } | |
| 223 | |||
| 224 | /***********************************************************************************/ | ||
| 225 | /***********************************************************************************/ | ||
| 226 | |||
| 227 | 42 | void NSModel::setExternalVec() | |
| 228 | { | ||
| 229 | // Retrieve instance | ||
| 230 | 42 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 231 | |||
| 232 | //========= | ||
| 233 | // ALE | ||
| 234 | //========= | ||
| 235 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 24 times.
|
42 | if ( r_instance.useALEformulation ) { |
| 236 | //1) Initialization of the vector containing the displacement of the mesh d^{n}, sequential std::vector | ||
| 237 | 18 | m_dispMesh.duplicateFrom(m_pSolverHE->sequentialSolution()); | |
| 238 | 18 | m_dispMesh.zeroEntries(); | |
| 239 | |||
| 240 | // fsi challenge | ||
| 241 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 12 times.
|
18 | if ( r_instance.useElasticExtension ) { |
| 242 | 6 | m_dispMesho_pp.duplicateFrom(m_pSolverHE->solution()); | |
| 243 | 6 | m_dispMesho_pp.zeroEntries(); | |
| 244 | } | ||
| 245 | |||
| 246 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
|
18 | if ( Tools::equal(r_instance.theta,0.5) ) { |
| 247 | ✗ | m_dispMesho.duplicateFrom(m_pSolverHE->sequentialSolution()); | |
| 248 | ✗ | m_dispMesho.zeroEntries(); | |
| 249 | } | ||
| 250 | |||
| 251 | //2) Initialization of the std::vector containing the displacement of the mesh d^{n-1}, parallel std::vector | ||
| 252 | // m_solDispl == vectors of dispalcement at every time step (n>1) | ||
| 253 | // m_solDispl_0 == vectors of dispalcement at first time step (n==1) | ||
| 254 | // Two vectors of displacements are necessary because when a | ||
| 255 | // cardiac cycle re-start we need to read the std::vector of | ||
| 256 | // displacement at time 0, which is stored in m_solDispl_0 | ||
| 257 | 18 | m_solDispl.duplicateFrom(m_pSolverHE->solution()); | |
| 258 | 18 | m_solDispl_0.duplicateFrom(m_pSolverHE->solution()); | |
| 259 | 18 | m_solDispl_0.zeroEntries(); | |
| 260 | 18 | m_solDispl.zeroEntries(); | |
| 261 | |||
| 262 | //3) Initialization of the vector containing the velocity of the mesh w, sequential vector | ||
| 263 | 18 | m_velMesh.duplicateFrom(m_pSolverHE->sequentialSolution()); | |
| 264 | 18 | m_velMesh.zeroEntries(); | |
| 265 | |||
| 266 | // Passing the ao and dof | ||
| 267 | 18 | m_pSolverNS->pushBackExternalAO(m_pSolverHE->ao()); | |
| 268 | 18 | m_pSolverNS->pushBackExternalDof(m_pSolverHE->dof()); | |
| 269 | // Passing the velocity of the mesh to the NS problem | ||
| 270 | 18 | m_pSolverNS->pushBackExternalVec(m_velMesh); //externalVec(0) in m_linearProblem[0] | |
| 271 | // Resizing the auxiliar vectors needed in updateMesh() | ||
| 272 | 18 | m_tmpDisp.resize(m_numDofExtensionProblem); | |
| 273 | |||
| 274 | 18 | m_displTimeStep = 0.; | |
| 275 | |||
| 276 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 14 times.
|
18 | if (r_instance.nonLinearFluid) { |
| 277 | // create a sequential PetscVector whose size is equal to the total number of fluid elements | ||
| 278 | 4 | m_solPrev.duplicateFrom(m_pSolverNS->solution()); | |
| 279 | 4 | m_solPrev.zeroEntries(); | |
| 280 | // push it back to the sets of linear problems | ||
| 281 | 4 | m_pSolverNS->pushBackExternalVec(m_solPrev); | |
| 282 | } | ||
| 283 | } | ||
| 284 | |||
| 285 |
4/4✓ Branch 0 taken 8 times.
✓ Branch 1 taken 34 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 4 times.
|
42 | if ( r_instance.hasImmersedData && r_instance.hasDivDivStab ) { |
| 286 | |||
| 287 | // create a sequential PetscVector whose size is equal to the total number of fluid elements | ||
| 288 | 4 | int ivar = m_pSolverNS->listVariable().getVariableIdList(velocity); | |
| 289 | 4 | int imesh = m_pSolverNS->listVariable()[ivar].idMesh(); | |
| 290 | |||
| 291 |
2/4✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
|
4 | m_delta.createSeq(PETSC_COMM_SELF, m_pSolverNS->mesh(imesh)->getNumDomainElement() ); |
| 292 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | m_delta.set(0.); |
| 293 | |||
| 294 | // push it back to the sets of linear problems | ||
| 295 | 4 | m_pSolverNS->pushBackExternalVec(m_delta); // TODO D.C. not really smart to do a pushBack in vector since the order may change depending on the problem we solve... mayba a map is better... | |
| 296 | } | ||
| 297 | |||
| 298 | // for monolithic fsi | ||
| 299 | // m_dispSeq is a sequential PetscVector whose size is equal to the total number of fluid elements | ||
| 300 | // m_dispSeq.createSeq(PETSC_COMM_SELF, m_pSolverNS->numDof() ); | ||
| 301 | // m_dispSeq.set(0.); | ||
| 302 | // m_pSolverNS->pushBackExternalVec(m_dispSeq); //externalVec(0) in m_pSolverNS or 1 if ALE | ||
| 303 | 42 | } | |
| 304 | |||
| 305 | /***********************************************************************************/ | ||
| 306 | /***********************************************************************************/ | ||
| 307 | |||
| 308 | 60 | void NSModel::preAssemblingMatrixRHS(std::size_t iProblem) | |
| 309 | { | ||
| 310 | // Retrieve instance | ||
| 311 | 60 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 312 | |||
| 313 |
2/2✓ Branch 0 taken 42 times.
✓ Branch 1 taken 18 times.
|
60 | if ( iProblem == 0 ) { |
| 314 | 42 | m_solExtrapolate.duplicateFrom(m_pSolverNS->vector()); | |
| 315 | 42 | m_solExtrapolate.zeroEntries(); | |
| 316 | |||
| 317 | 42 | m_pSolverNS->initExtrapol(m_solExtrapolate); | |
| 318 | |||
| 319 | 42 | m_pSolverNS->solution().zeroEntries(); | |
| 320 | |||
| 321 | //! First assembly loop in iteration 0 to build static matrix. | ||
| 322 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 18 times.
|
42 | if ( r_instance.useALEformulation == 0 ) { |
| 323 | |||
| 324 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
|
24 | if ( r_instance.NSStabType == 1 ) |
| 325 | ✗ | m_pSolverNS->assembleFaceOrientedStabilization(); | |
| 326 | |||
| 327 | 24 | m_pSolverNS->assembleMatrixRHS(MpiInfo::rankProc()); | |
| 328 | } | ||
| 329 | } | ||
| 330 | |||
| 331 | //========= | ||
| 332 | // ALE | ||
| 333 | //========= | ||
| 334 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 24 times.
|
60 | if ( r_instance.useALEformulation ) { |
| 335 |
4/4✓ Branch 0 taken 18 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 6 times.
|
36 | if ( (iProblem == 1) && (!r_instance.useElasticExtension) ) { |
| 336 | // First assembly loop in iteration 0 to build static matrix of the extension problem | ||
| 337 | 12 | m_pSolverHE->assembleMatrixRHS(MpiInfo::rankProc()); | |
| 338 | // Save the static part of the extension problem in the matrix _A | ||
| 339 | 12 | m_pSolverHE->copyMatrixRHS(); | |
| 340 | } | ||
| 341 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
36 | if ( (r_instance.useALEformulation == 2) && (iProblem == 0) ) { |
| 342 | // assembly loop to put values at the locations of the complementary pattern | ||
| 343 | ✗ | if ( r_instance.NSStabType == 1 ) | |
| 344 | ✗ | m_pSolverNS->assembleFaceOrientedStabilization(); | |
| 345 | |||
| 346 | // First assembly loop in iteration 0 to build the initial mass matrix (initial configuration) | ||
| 347 | ✗ | m_pSolverNS->assembleMatrixRHS(MpiInfo::rankProc()); | |
| 348 | } | ||
| 349 | } | ||
| 350 | 60 | } | |
| 351 | |||
| 352 | /***********************************************************************************/ | ||
| 353 | /***********************************************************************************/ | ||
| 354 | |||
| 355 | 1904 | void NSModel::postAssemblingMatrixRHS(std::size_t iProblem) | |
| 356 | { | ||
| 357 | // Retrieve instance | ||
| 358 | 1904 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 359 | |||
| 360 | // matrix(0) = matrix(0) + matrix(1) (add static matrix to the dynamic matrix to build | ||
| 361 | // complete matrix of the system ) | ||
| 362 | 1904 | m_linearProblem[iProblem]->addMatrixRHS(); | |
| 363 | |||
| 364 |
3/4✓ Branch 0 taken 1544 times.
✓ Branch 1 taken 360 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1544 times.
|
1904 | if ( iProblem == 0 && r_instance.withCVG ) |
| 365 | ✗ | m_pSolverNS->createAndCopyMatrixRHSWithoutBC(); | |
| 366 | |||
| 367 | //Apply CVGraphInterface condition directly in RHS | ||
| 368 |
3/4✓ Branch 0 taken 1544 times.
✓ Branch 1 taken 360 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1544 times.
|
1904 | if ( iProblem == 0 && r_instance.CVGraphInterfaceVariationalBC == 1 ) { |
| 369 | ✗ | m_pSolverNS->vector().setValues(m_cvgraphInterf->SizeListOfDofForVariable(), | |
| 370 | ✗ | m_cvgraphInterf->listOfDofInInterface(), m_cvgraphInterf->ComputedVariableToImpose(), ADD_VALUES); | |
| 371 | } | ||
| 372 | 1904 | } | |
| 373 | |||
| 374 | /***********************************************************************************/ | ||
| 375 | /***********************************************************************************/ | ||
| 376 | |||
| 377 | ✗ | void NSModel::RISforward() | |
| 378 | { | ||
| 379 | // Retrieve instance | ||
| 380 | ✗ | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 381 | |||
| 382 | ✗ | m_chek_use_RISforward = true; | |
| 383 | |||
| 384 | ✗ | if ( m_risModel->admissibleStatus() == 1 ) { | |
| 385 | ✗ | writeSolution(); | |
| 386 | ✗ | if ( m_lumpedModelBC && MpiInfo::rankProc() == 0 ) | |
| 387 | ✗ | m_lumpedModelBC->write(r_instance.resultDir+"/output_windkessel"); | |
| 388 | |||
| 389 | ✗ | if ( MpiInfo::rankProc() == 0 ) | |
| 390 | ✗ | m_risModel->write(r_instance.resultDir+"/outputValves"); | |
| 391 | |||
| 392 | // Advance time step. | ||
| 393 | ✗ | updateTime(); | |
| 394 | |||
| 395 | ✗ | if ( r_instance.useALEformulation == 2 ) { | |
| 396 | |||
| 397 | ✗ | if ( m_fstransient->iteration == 1 ) { | |
| 398 | //create m_storedMatrix with the same structure of matrix(1) | ||
| 399 | ✗ | m_storedMatrix.duplicateFrom(m_pSolverNS->matrix(1), MAT_COPY_VALUES); | |
| 400 | } | ||
| 401 | |||
| 402 | ✗ | if ( m_fstransient->iteration > 1 ) { | |
| 403 | // add Matrix(1) to matrix(0) | ||
| 404 | ✗ | m_pSolverNS->matrix(0).axpy(1,m_pSolverNS->matrix(1),SAME_NONZERO_PATTERN); | |
| 405 | // copy matrix(1) in m_storedMatrix to use it in the NS recomputation | ||
| 406 | ✗ | m_storedMatrix.copyFrom(m_pSolverNS->matrix(1), SAME_NONZERO_PATTERN); | |
| 407 | //clear matrix(1) to be updated in the forward() | ||
| 408 | ✗ | m_pSolverNS->clearMatrix(1); | |
| 409 | } | ||
| 410 | } | ||
| 411 | |||
| 412 | // This is usefeul to read data_displacements | ||
| 413 | ✗ | m_displTimeStep += m_fstransient->timeStep; | |
| 414 | // Print time information | ||
| 415 | ✗ | printNewTimeIterationBanner(); | |
| 416 | |||
| 417 | ✗ | m_risModel->UpdateResistances(); | |
| 418 | |||
| 419 | ✗ | if(r_instance.CardiacCycle) | |
| 420 | ✗ | m_cardiacCycle->ApplyCardiacCycleInput(); | |
| 421 | |||
| 422 | ✗ | forward(); | |
| 423 | ✗ | m_risModel->CheckValveStatus(); | |
| 424 | } | ||
| 425 | |||
| 426 | ✗ | while( m_risModel->admissibleStatus() == 0 ) { // This "while" is important, an "if" does not work in all cases | |
| 427 | ✗ | PetscPrintf(MpiInfo::petscComm(),"\n Recomputation with the new RIS status\n"); | |
| 428 | |||
| 429 | ✗ | for(std::size_t ilr=0; ilr<m_linearProblem.size(); ilr++) { | |
| 430 | ✗ | m_linearProblem[ilr]->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs); // 0 is flagMatrixRHS clear matrix(0) and RHS | |
| 431 | } | ||
| 432 | |||
| 433 | ✗ | if( (r_instance.useALEformulation == 2) && ( m_fstransient->iteration > 1 ) ) { | |
| 434 | // we recomputed the matrix(0) at t^{n} and we add m_storedMatrix == matrix(1) a t^{n-1} | ||
| 435 | ✗ | m_pSolverNS->matrix(0).axpy(1,m_storedMatrix,SAME_NONZERO_PATTERN); | |
| 436 | ✗ | m_storedMatrix.setUnfactored(); | |
| 437 | ✗ | m_storedMatrix.zeroEntries(); | |
| 438 | } | ||
| 439 | |||
| 440 | ✗ | m_risModel->UpdateResistances(); | |
| 441 | |||
| 442 | ✗ | if( r_instance.CardiacCycle ) | |
| 443 | ✗ | m_cardiacCycle->ApplyCardiacCycleInput(); | |
| 444 | |||
| 445 | ✗ | forward(); | |
| 446 | ✗ | m_risModel->CheckValveStatus(); | |
| 447 | } | ||
| 448 | } | ||
| 449 | |||
| 450 | /***********************************************************************************/ | ||
| 451 | /***********************************************************************************/ | ||
| 452 | |||
| 453 | 1444 | void NSModel::forward() | |
| 454 | { | ||
| 455 | 1444 | prepareForward(); | |
| 456 | 1444 | solveNS(); | |
| 457 | 1444 | } | |
| 458 | |||
| 459 | /***********************************************************************************/ | ||
| 460 | /***********************************************************************************/ | ||
| 461 | |||
| 462 | 1544 | void NSModel::prepareForward() | |
| 463 | { | ||
| 464 | // Retrieve instance | ||
| 465 | 1544 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 466 | |||
| 467 | // Initial displacement for ALE | ||
| 468 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 1544 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1544 times.
|
1544 | if ( r_instance.initALEDispByFile && (m_fstransient->iteration == 0 ) ) { |
| 469 | |||
| 470 | ✗ | m_pSolverHE->readMatchFile_ALE(m_pSolverHE->listVariable(), r_instance.MoveMeshMatchFile); | |
| 471 | ✗ | m_pSolverHE->identifyIdDofToMatch_ALE(m_pSolverHE->dof()); | |
| 472 | ✗ | m_pSolverHE->readDataDisplacement(m_ios, 0.0); | |
| 473 | int low, high; | ||
| 474 | ✗ | VecGetOwnershipRange(m_pSolverHE->solution().toPetsc(),&low,&high); | |
| 475 | ✗ | int size = high-low; | |
| 476 | ✗ | std::vector<double> vv(size); | |
| 477 | ✗ | std::vector<int> ix(size); | |
| 478 | ✗ | std::iota(ix.begin(),ix.end(),low); // fills ix with values from low to high | |
| 479 | |||
| 480 | ✗ | VecGetValues(m_pSolverHE->HarmExtSol().toPetsc(),size,ix.data(),vv.data()); | |
| 481 | |||
| 482 | ✗ | VecSetValues(m_pSolverHE->solution().toPetsc(),size,ix.data(),vv.data(), INSERT_VALUES); | |
| 483 | ✗ | m_pSolverHE->solution().assembly(); | |
| 484 | ✗ | m_pSolverHE->gatherSolution(); | |
| 485 | ✗ | updateMesh(); | |
| 486 | } | ||
| 487 | |||
| 488 |
1/2✓ Branch 0 taken 1544 times.
✗ Branch 1 not taken.
|
1544 | if ( !r_instance.RISModels ) { |
| 489 | |||
| 490 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1544 times.
|
1544 | if ( r_instance.useFDlagrangeMult ) { |
| 491 | |||
| 492 | // Compute the displacements of the interface meshes | ||
| 493 | ✗ | m_pSolverNS->computeItfDisp(*m_itfDisp); | |
| 494 | |||
| 495 | // Move the interface meshes | ||
| 496 | ✗ | m_pSolverNS->moveItfMeshes(); | |
| 497 | |||
| 498 | // Write displacement and normal | ||
| 499 | ✗ | const int ivar = m_pSolverNS->listVariable().getVariableIdList(lagMultiplier); | |
| 500 | ✗ | const int imesh = m_pSolverNS->listVariable()[ivar].idMesh(); | |
| 501 | ✗ | writeCurrentDisplacement(imesh, "displacement"); | |
| 502 | ✗ | writeCurrentNormalVector(); | |
| 503 | } | ||
| 504 | |||
| 505 | // Write solution for postprocessing (if required) | ||
| 506 | 1544 | writeSolution(); | |
| 507 |
2/2✓ Branch 0 taken 40 times.
✓ Branch 1 taken 1504 times.
|
1544 | if ( m_lumpedModelBC ) { |
| 508 |
2/2✓ Branch 1 taken 10 times.
✓ Branch 2 taken 30 times.
|
40 | if ( MpiInfo::rankProc() == 0 ) |
| 509 |
1/2✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
|
10 | m_lumpedModelBC->write(r_instance.resultDir+"/output_windkessel"); |
| 510 | } | ||
| 511 | |||
| 512 | 1544 | updateTime(); //Here the matrix(0), the non constant one, is cleared and the rhs also. | |
| 513 | |||
| 514 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 1544 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1544 times.
|
1544 | if( (r_instance.useALEformulation == 2) && ( m_fstransient->iteration > 1 ) ) { |
| 515 | ✗ | m_pSolverNS->matrix(0).axpy(1,m_pSolverNS->matrix(1),SAME_NONZERO_PATTERN); | |
| 516 | ✗ | m_pSolverNS->clearMatrix(1); | |
| 517 | } | ||
| 518 | |||
| 519 | // Print time information | ||
| 520 | 1544 | printNewTimeIterationBanner(); | |
| 521 | |||
| 522 | // This is usefeul to read data_displacements | ||
| 523 | 1544 | m_displTimeStep += m_fstransient->timeStep; | |
| 524 | } | ||
| 525 | ✗ | else if ( /*r_instance.RISModels &&*/ m_chek_use_RISforward == false ) { | |
| 526 | |||
| 527 | ✗ | FEL_ERROR("Your are not using the corret forward for valve implementation !") | |
| 528 | } | ||
| 529 | |||
| 530 | 1544 | m_pSolverNS->onPreviousMesh(); | |
| 531 | |||
| 532 | // Function only at iteration == 1 | ||
| 533 | 1544 | firstIteration(); | |
| 534 | |||
| 535 |
2/2✓ Branch 1 taken 1502 times.
✓ Branch 2 taken 42 times.
|
1544 | if ( m_fstransient->iteration > 1 ) { |
| 536 | 1502 | m_bdfNS.update(m_pSolverNS->solution()); | |
| 537 | 1502 | m_bdfNS.extrapolate(m_solExtrapolate); | |
| 538 | 1502 | m_pSolverNS->updateExtrapol(m_solExtrapolate); | |
| 539 | |||
| 540 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1502 times.
|
1502 | if ( Tools::equal(r_instance.theta,0.5) ) { |
| 541 | ✗ | m_CrankNicolsonNS->update(m_pSolverNS->solution()); | |
| 542 | ✗ | m_CrankNicolsonNS->extrapolate(m_solExtrapolate); | |
| 543 | ✗ | m_pSolverNS->updateExtrapol(m_solExtrapolate); | |
| 544 | } | ||
| 545 | } | ||
| 546 | |||
| 547 | // If lumpedModel bc computation | ||
| 548 | 1544 | computeLumpedModelBC(); | |
| 549 | |||
| 550 | 1544 | m_bdfNS.computeRHSTime(m_fstransient->timeStep); | |
| 551 | 1544 | m_pSolverNS->gatherVectorBeforeAssembleMatrixRHS(); | |
| 552 | |||
| 553 | // If fictitious domain + penalization update interface mesh | ||
| 554 |
2/2✓ Branch 0 taken 320 times.
✓ Branch 1 taken 1224 times.
|
1544 | if ( r_instance.hasImmersedData ) { |
| 555 | 320 | GeometricMeshRegion& mesh = *m_pSolverNS->mesh(); | |
| 556 | 320 | m_immersed.updateInterface(mesh); // TODO D.C. make more general | |
| 557 | |||
| 558 |
2/2✓ Branch 0 taken 200 times.
✓ Branch 1 taken 120 times.
|
320 | if ( r_instance.hasDivDivStab ) { |
| 559 | // call updateIntersection of the model "Immersed" in order to set the values of the PetscVector m_delta for the main process and to broadcast it to the other processes | ||
| 560 | 200 | m_immersed.updateIntersection(mesh, m_delta, MpiInfo::petscComm()); | |
| 561 | } | ||
| 562 | } | ||
| 563 | |||
| 564 | // If fictitious domain with lagrange multiplier | ||
| 565 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1544 times.
|
1544 | if ( r_instance.useFDlagrangeMult ) { |
| 566 | |||
| 567 | // Localize interface quadrature point in fluid mesh and evaluate new pattern | ||
| 568 | ✗ | m_pSolverNS->computeSolidFluidMaps(); | |
| 569 | ✗ | m_pSolverNS->evaluateFluidStructurePattern(); | |
| 570 | |||
| 571 | ✗ | m_pSolverNS->computeLagrangeMultiplierOnInterface(FlagMatrixRHS::only_matrix); | |
| 572 | |||
| 573 | ✗ | if ( r_instance.BHstab ) | |
| 574 | ✗ | m_pSolverNS->assembleBarbosaHughesStab(); | |
| 575 | |||
| 576 | ✗ | if ( r_instance.BPstab ) | |
| 577 | ✗ | m_pSolverNS->assembleBrezziPitkarantaStab(); | |
| 578 | |||
| 579 | ✗ | if ( r_instance.massConstraint ) { | |
| 580 | ✗ | m_pSolverNS->assembleJapanConstraintBoundary(); | |
| 581 | ✗ | m_pSolverNS->assembleJapanConstraintInterface(); | |
| 582 | } | ||
| 583 | } | ||
| 584 | |||
| 585 |
2/2✓ Branch 0 taken 360 times.
✓ Branch 1 taken 1184 times.
|
1544 | if ( r_instance.useALEformulation ) { |
| 586 | |||
| 587 | 360 | computeALEterms(); | |
| 588 | } else { | ||
| 589 | |||
| 590 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1184 times.
|
1184 | if ( r_instance.NSStabType == 1 ) |
| 591 | ✗ | m_pSolverNS->assembleFaceOrientedStabilization(); | |
| 592 | |||
| 593 | 1184 | m_pSolverNS->assembleMatrixRHS(MpiInfo::rankProc()); | |
| 594 | } | ||
| 595 | |||
| 596 | // TODO: we can use the functions createAndCopyMatrixRHSWithoutBC and computeResidual in LinearProblem | ||
| 597 |
2/2✓ Branch 0 taken 100 times.
✓ Branch 1 taken 1444 times.
|
1544 | if ( r_instance.fsiCoupling ) |
| 598 | 100 | m_rhsWihoutBC.copyFrom(m_pSolverNS->vector()); | |
| 599 | |||
| 600 | // Specific operations before solve the system. | ||
| 601 | 1544 | postAssemblingMatrixRHS(0); | |
| 602 | |||
| 603 | // TODO: we can use the functions createAndCopyMatrixRHSWithoutBC and computeResidual in LinearProblem | ||
| 604 |
3/4✓ Branch 0 taken 100 times.
✓ Branch 1 taken 1444 times.
✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
|
1544 | if ( r_instance.fsiCoupling && !r_instance.useFDlagrangeMult ) { |
| 605 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 94 times.
|
100 | if ( m_matrixWithoutBCcreated == 0 ) { |
| 606 | 6 | m_matrixWithoutBC.duplicateFrom(m_pSolverNS->matrix(),MAT_COPY_VALUES); | |
| 607 | 6 | m_matrixWithoutBCcreated = true; | |
| 608 | } | ||
| 609 | 100 | m_matrixWithoutBC.copyFrom(m_pSolverNS->matrix(), SAME_NONZERO_PATTERN); | |
| 610 | } | ||
| 611 | 1544 | } | |
| 612 | |||
| 613 | /***********************************************************************************/ | ||
| 614 | /***********************************************************************************/ | ||
| 615 | |||
| 616 | 1544 | void NSModel::solveNS() | |
| 617 | { | ||
| 618 | 1544 | const auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 619 | |||
| 620 | 1544 | const int ivar = m_pSolverNS->listVariable().getVariableIdList(velocity); | |
| 621 | 1544 | const int imesh = m_pSolverNS->listVariable()[ivar].idMesh(); | |
| 622 | |||
| 623 |
4/6✓ Branch 0 taken 100 times.
✓ Branch 1 taken 1444 times.
✓ Branch 2 taken 100 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
|
1544 | if ( r_instance.fsiCoupling && !r_instance.hasImmersedData && !r_instance.useFDlagrangeMult ) { |
| 624 |
2/2✓ Branch 1 taken 6 times.
✓ Branch 2 taken 94 times.
|
100 | if ( m_fstransient->iteration == 1 ) { |
| 625 |
1/2✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
|
6 | if ( m_hasLumpedMassRN ) { |
| 626 | 6 | PetscScalar dti = 1./m_fstransient->timeStep; | |
| 627 | 6 | m_lumpedMassRN.scale(dti); | |
| 628 | } | ||
| 629 | } | ||
| 630 | |||
| 631 |
3/4✓ Branch 0 taken 70 times.
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70 times.
|
100 | if ( r_instance.bcInRefConfiguration && !m_hasLumpedMassRN ) { |
| 632 | ✗ | m_pSolverNS->meshLocal(imesh)->moveMesh(m_dispMeshVector,0.0); | |
| 633 | } | ||
| 634 |
3/4✓ Branch 0 taken 70 times.
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70 times.
|
100 | if ( r_instance.bcInRefConfiguration && r_instance.readPreStressedDisplFromFile ) { |
| 635 | ✗ | m_pSolverNS->meshLocal(imesh)->moveMesh(m_preStressedDispMeshVector,1.0); | |
| 636 | } | ||
| 637 |
3/6✓ Branch 0 taken 260 times.
✓ Branch 1 taken 1184 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 260 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
1544 | } else if ( r_instance.useALEformulation && r_instance.bcInRefConfiguration && r_instance.NSequationFlag != 0 ) { |
| 638 | ✗ | m_pSolverNS->meshLocal(imesh)->moveMesh(m_dispMeshVector,0.0); | |
| 639 | } | ||
| 640 | |||
| 641 |
2/2✓ Branch 1 taken 300 times.
✓ Branch 2 taken 1244 times.
|
1544 | if ( FelisceParam::verbose() ) { |
| 642 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 300 times.
|
300 | if ( r_instance.nonLinearFluid ) |
| 643 | ✗ | PetscPrintf(MpiInfo::petscComm(),"NS non-linear problem \n"); | |
| 644 | else | ||
| 645 | 300 | PetscPrintf(MpiInfo::petscComm(),"NS linear problem \n"); | |
| 646 | } | ||
| 647 | |||
| 648 |
2/2✓ Branch 0 taken 320 times.
✓ Branch 1 taken 1224 times.
|
1544 | if ( r_instance.hasImmersedData ) { |
| 649 | 320 | AO& ao = m_pSolverNS->ao(); | |
| 650 | 320 | PetscMatrix& matrix = m_pSolverNS->matrix(0); | |
| 651 | 320 | PetscVector& rhs = m_pSolverNS->vector(); | |
| 652 | 320 | Dof& dof = m_pSolverNS->dof(); | |
| 653 | 320 | int idVarVel = m_pSolverNS->idVarVel(); | |
| 654 | 320 | m_immersed.applyPenaltyToMatrixAndRHS(ao,dof,idVarVel,matrix,rhs); | |
| 655 | } | ||
| 656 | |||
| 657 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1544 times.
|
1544 | if ( r_instance.useFDlagrangeMult) |
| 658 | ✗ | m_pSolverNS->computeLagrangeMultiplierOnInterface(FlagMatrixRHS::only_rhs); | |
| 659 | |||
| 660 |
2/2✓ Branch 1 taken 300 times.
✓ Branch 2 taken 1244 times.
|
1544 | if ( FelisceParam::verbose() ) |
| 661 | 300 | PetscPrintf(MpiInfo::petscComm(),"Applying BC \n"); | |
| 662 | |||
| 663 | //Apply essential transient boundary conditions. | ||
| 664 | 1544 | m_pSolverNS->finalizeEssBCTransient(); | |
| 665 | |||
| 666 | #ifdef FELISCE_WITH_CVGRAPH | ||
| 667 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 1544 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1544 times.
✗ Branch 7 not taken.
|
1544 | if ( !r_instance.withCVG || ! m_pSolverNS->slave()->redoTimeStepOnlyCVGData() ) { |
| 668 | 1544 | m_pSolverNS->applyBC(r_instance.essentialBoundaryConditionsMethod, MpiInfo::rankProc()); | |
| 669 | } else { | ||
| 670 | ✗ | m_pSolverNS->applyOnlyCVGBoundaryCondition(); | |
| 671 | } | ||
| 672 | #else | ||
| 673 | m_pSolverNS->applyBC(r_instance.essentialBoundaryConditionsMethod, MpiInfo::rankProc()); | ||
| 674 | #endif | ||
| 675 | |||
| 676 |
2/2✓ Branch 0 taken 100 times.
✓ Branch 1 taken 1444 times.
|
1544 | if ( r_instance.fsiCoupling ) { |
| 677 | |||
| 678 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 100 times.
|
100 | if ( r_instance.useLumpedMassRN ) { |
| 679 | ✗ | m_pSolverNS->matrix().diagonalSet(m_lumpedMassRN,ADD_VALUES); | |
| 680 | } | ||
| 681 |
3/4✓ Branch 0 taken 70 times.
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70 times.
|
100 | if ( r_instance.bcInRefConfiguration && r_instance.readPreStressedDisplFromFile ) { |
| 682 | ✗ | m_pSolverNS->meshLocal(imesh)->moveMesh(m_preStressedDispMeshVector,-1.0); | |
| 683 | } | ||
| 684 |
3/4✓ Branch 0 taken 70 times.
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70 times.
|
100 | if ( r_instance.bcInRefConfiguration && !m_hasLumpedMassRN ) { |
| 685 | ✗ | m_pSolverNS->meshLocal(imesh)->moveMesh(m_dispMeshVector,1.0); | |
| 686 | } | ||
| 687 | |||
| 688 | 100 | m_pSolverNS->kspInterface().setDiagonalScale(); | |
| 689 | } | ||
| 690 |
3/6✓ Branch 0 taken 260 times.
✓ Branch 1 taken 1184 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 260 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
1444 | else if ( r_instance.useALEformulation && r_instance.bcInRefConfiguration && r_instance.NSequationFlag != 0 ) { |
| 691 | ✗ | m_pSolverNS->meshLocal(imesh)->moveMesh(m_dispMeshVector,1.0); | |
| 692 | } | ||
| 693 | |||
| 694 |
2/2✓ Branch 1 taken 300 times.
✓ Branch 2 taken 1244 times.
|
1544 | if ( FelisceParam::verbose() ) |
| 695 | 300 | PetscPrintf(MpiInfo::petscComm(),"Solving NS system \n"); | |
| 696 | |||
| 697 |
2/2✓ Branch 0 taken 960 times.
✓ Branch 1 taken 584 times.
|
1544 | if ( r_instance.nonLinearFluid ) { |
| 698 | // solve non Linear system | ||
| 699 | 960 | m_pSolverNS->iterativeSolve(MpiInfo::rankProc(), MpiInfo::numProc(), ApplyNaturalBoundaryConditions::no, FlagMatrixRHS::matrix_and_rhs); | |
| 700 | } else { | ||
| 701 | // solve Linear system | ||
| 702 | 584 | m_pSolverNS->solve(MpiInfo::rankProc(), MpiInfo::numProc()); | |
| 703 | } | ||
| 704 | |||
| 705 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1544 times.
|
1544 | if ( r_instance.withCVG ) { |
| 706 | ✗ | m_pSolverNS->computeResidual(); | |
| 707 | } | ||
| 708 | |||
| 709 |
2/2✓ Branch 0 taken 100 times.
✓ Branch 1 taken 1444 times.
|
1544 | if ( r_instance.fsiCoupling ) { |
| 710 | 100 | m_pSolverNS->gatherSolution(); | |
| 711 | |||
| 712 |
1/2✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
|
100 | if ( FelisceParam::verbose() ) |
| 713 | 100 | PetscPrintf(MpiInfo::petscComm(),"Evaluating residual \n"); | |
| 714 | |||
| 715 | // TODO: we can use the functions createAndCopyMatrixRHSWithoutBC and computeResidual in LinearProblemNS | ||
| 716 | // Compute residual | ||
| 717 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 100 times.
|
100 | if ( r_instance.useFDlagrangeMult ) { |
| 718 | ✗ | m_pSolverNS->computeInterfaceStress(); | |
| 719 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 100 times.
|
100 | } else if ( r_instance.hasImmersedData ) { |
| 720 | ✗ | AO& ao = m_pSolverNS->ao(); | |
| 721 | ✗ | PetscVector& seqVel = m_pSolverNS->sequentialSolution(); | |
| 722 | ✗ | Dof& dof = m_pSolverNS->dof(); | |
| 723 | ✗ | int idVarVel = m_pSolverNS->idVarVel(); | |
| 724 | ✗ | m_immersed.computePenaltyForce(ao,dof,idVarVel,seqVel); | |
| 725 | } else { | ||
| 726 | 100 | m_matrixWithoutBC.scale(-1.); | |
| 727 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 100 times.
|
100 | if ( r_instance.nonLinearFluid ) |
| 728 | ✗ | m_rhsWihoutBC.scale(-1.); | |
| 729 | |||
| 730 | 100 | multAdd(m_matrixWithoutBC,m_pSolverNS->solution(),m_rhsWihoutBC,m_residual); | |
| 731 | |||
| 732 | 100 | m_pSolverNS->gatherVector(m_residual, m_seqResidual); | |
| 733 | 100 | m_matrixWithoutBC.scale(-1.); | |
| 734 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 100 times.
|
100 | if ( r_instance.nonLinearFluid ) |
| 735 | ✗ | m_rhsWihoutBC.scale(-1.); | |
| 736 | } | ||
| 737 | } | ||
| 738 | |||
| 739 | // (interface) displacement update, monolthic FSI | ||
| 740 | // d^n = d^{n-1} + \tau * u^n | ||
| 741 | // m_disp.axpy(m_fstransient->timeStep,m_pSolverNS->solution()); | ||
| 742 | // m_pSolverNS->gatherVector(m_disp, m_dispSeq); | ||
| 743 | // m_pSolverNS->gatherSolution(); | ||
| 744 | 1544 | } | |
| 745 | |||
| 746 | /***********************************************************************************/ | ||
| 747 | /***********************************************************************************/ | ||
| 748 | |||
| 749 | 984 | void NSModel::solveLinearizedNS(bool linearFlag) | |
| 750 | { | ||
| 751 | // Retrieve instance | ||
| 752 | 984 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 753 | |||
| 754 | 984 | m_pSolverNS->linearizedFlag() = linearFlag; | |
| 755 | |||
| 756 | // back to original matrix | ||
| 757 | 984 | m_pSolverNS->kspInterface().setDiagonalScaleFix(); | |
| 758 | |||
| 759 |
2/2✓ Branch 0 taken 874 times.
✓ Branch 1 taken 110 times.
|
984 | if ( linearFlag ) |
| 760 | 874 | m_pSolverNS->vector().zeroEntries(); | |
| 761 | else | ||
| 762 | 110 | m_pSolverNS->vector().copyFrom(m_rhsWihoutBC); | |
| 763 | |||
| 764 |
5/6✓ Branch 0 taken 874 times.
✓ Branch 1 taken 110 times.
✓ Branch 3 taken 874 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 874 times.
✓ Branch 6 taken 110 times.
|
984 | if ( linearFlag && FelisceParam::verbose() ) |
| 765 | 874 | PetscPrintf(MpiInfo::petscComm(),"NS linearized solver \n"); | |
| 766 |
1/2✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
|
110 | else if ( FelisceParam::verbose() ) |
| 767 | 110 | PetscPrintf(MpiInfo::petscComm(),"NS solver \n"); | |
| 768 | |||
| 769 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 984 times.
|
984 | if ( r_instance.hasImmersedData ) { |
| 770 | ✗ | AO& ao = m_pSolverNS->ao(); | |
| 771 | ✗ | PetscMatrix& matrix = m_pSolverNS->matrix(0); | |
| 772 | ✗ | PetscVector& rhs = m_pSolverNS->vector(); | |
| 773 | ✗ | Dof& dof = m_pSolverNS->dof(); | |
| 774 | ✗ | int idVarVel = m_pSolverNS->idVarVel(); | |
| 775 | ✗ | m_immersed.applyPenaltyToMatrixAndRHS(ao,dof,idVarVel,matrix,rhs,true); | |
| 776 | } | ||
| 777 | |||
| 778 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 984 times.
|
984 | if ( r_instance.useFDlagrangeMult) { |
| 779 | ✗ | m_pSolverNS->computeLagrangeMultiplierOnInterface(FlagMatrixRHS::only_rhs); | |
| 780 | } | ||
| 781 | |||
| 782 |
1/2✓ Branch 1 taken 984 times.
✗ Branch 2 not taken.
|
984 | if (FelisceParam::verbose()) |
| 783 | 984 | PetscPrintf(MpiInfo::petscComm(),"Applying BC \n"); | |
| 784 | |||
| 785 | // Apply essential transient boundary conditions. | ||
| 786 | 984 | m_pSolverNS->finalizeEssBCTransient(); | |
| 787 | |||
| 788 |
2/2✓ Branch 0 taken 874 times.
✓ Branch 1 taken 110 times.
|
984 | if ( linearFlag ) // only Dirichlet BC on the interface is applied (RHS) |
| 789 | 874 | m_pSolverNS->applyBC(r_instance.essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs, 0, true, ApplyNaturalBoundaryConditions::no); | |
| 790 | else // all BC applied (RHS) | ||
| 791 | 110 | m_pSolverNS->applyBC(r_instance.essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs, 0, true, ApplyNaturalBoundaryConditions::yes); | |
| 792 | |||
| 793 | 984 | m_pSolverNS->kspInterface().setDiagonalScale(); | |
| 794 | |||
| 795 |
1/2✓ Branch 1 taken 984 times.
✗ Branch 2 not taken.
|
984 | if ( FelisceParam::verbose() ) { |
| 796 |
2/2✓ Branch 0 taken 874 times.
✓ Branch 1 taken 110 times.
|
984 | if ( linearFlag ) |
| 797 | 874 | PetscPrintf(MpiInfo::petscComm(),"Solving linearized NS system (same matrix) \n"); | |
| 798 | else | ||
| 799 | 110 | PetscPrintf(MpiInfo::petscComm(),"Solving NS system (same matrix) \n"); | |
| 800 | } | ||
| 801 | |||
| 802 | // Solve linear system with same matrix as in solveNS | ||
| 803 | 984 | m_pSolverNS->solveWithSameMatrix(); | |
| 804 | |||
| 805 | 984 | m_pSolverNS->gatherSolution(); | |
| 806 | |||
| 807 |
1/2✓ Branch 1 taken 984 times.
✗ Branch 2 not taken.
|
984 | if ( FelisceParam::verbose() ) |
| 808 | 984 | PetscPrintf(MpiInfo::petscComm(),"Evaluating residual \n"); | |
| 809 | |||
| 810 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 984 times.
|
984 | if ( r_instance.useFDlagrangeMult ) { |
| 811 | ✗ | m_pSolverNS->computeInterfaceStress(); | |
| 812 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 984 times.
|
984 | } else if ( r_instance.hasImmersedData ) { |
| 813 | ✗ | AO& ao = m_pSolverNS->ao(); | |
| 814 | ✗ | PetscVector& seqVel = m_pSolverNS->sequentialSolution(); | |
| 815 | ✗ | Dof& dof = m_pSolverNS->dof(); | |
| 816 | ✗ | int idVarVel = m_pSolverNS->idVarVel(); | |
| 817 | ✗ | m_immersed.computePenaltyForce(ao,dof,idVarVel,seqVel); | |
| 818 | } else { | ||
| 819 | //Compute residual | ||
| 820 | 984 | m_matrixWithoutBC.scale(-1.); | |
| 821 |
2/2✓ Branch 0 taken 874 times.
✓ Branch 1 taken 110 times.
|
984 | if ( linearFlag ) { |
| 822 | 874 | m_pSolverNS->vector().zeroEntries(); | |
| 823 | 874 | multAdd(m_matrixWithoutBC,m_pSolverNS->solution(),m_pSolverNS->vector(),m_residual); | |
| 824 | } | ||
| 825 | else | ||
| 826 | 110 | multAdd(m_matrixWithoutBC,m_pSolverNS->solution(),m_rhsWihoutBC,m_residual); | |
| 827 | 984 | m_pSolverNS->gatherVector(m_residual, m_seqResidual); | |
| 828 | 984 | m_matrixWithoutBC.scale(-1.); | |
| 829 | } | ||
| 830 | |||
| 831 | 984 | m_pSolverNS->linearizedFlag() = false; | |
| 832 | 984 | } | |
| 833 | |||
| 834 | /***********************************************************************************/ | ||
| 835 | /***********************************************************************************/ | ||
| 836 | |||
| 837 | 1544 | void NSModel::computeLumpedModelBC() | |
| 838 | { | ||
| 839 | // Retrieve instance | ||
| 840 | 1544 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 841 | |||
| 842 |
2/2✓ Branch 0 taken 40 times.
✓ Branch 1 taken 1504 times.
|
1544 | if (m_lumpedModelBC) { |
| 843 | // update the fluxes | ||
| 844 | 40 | m_pSolverNS->computeMeanQuantity(velocity, m_labelFluxListLumpedModelBC, m_fluxLumpedModelBC); | |
| 845 | |||
| 846 |
2/2✓ Branch 1 taken 80 times.
✓ Branch 2 taken 40 times.
|
120 | for(std::size_t i=0; i<m_lumpedModelBC->size(); i++) { |
| 847 | 80 | m_lumpedModelBC->Q(i) = m_fluxLumpedModelBC[i] ; | |
| 848 | } | ||
| 849 | |||
| 850 | // update compliance if non linear compliance | ||
| 851 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
|
40 | if ( r_instance.lumpedModelBC_C_type == 2 ) { |
| 852 | ✗ | m_pSolverNS->updateVolume(); | |
| 853 | } | ||
| 854 | |||
| 855 | // iterate the lumpedModelBC parameters | ||
| 856 | 40 | m_lumpedModelBC->iterate(); | |
| 857 | } | ||
| 858 | 1544 | } | |
| 859 | |||
| 860 | /***********************************************************************************/ | ||
| 861 | /***********************************************************************************/ | ||
| 862 | |||
| 863 | 1544 | void NSModel::firstIteration() | |
| 864 | { | ||
| 865 | // Retrieve instance | ||
| 866 | 1544 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 867 | |||
| 868 |
2/2✓ Branch 1 taken 42 times.
✓ Branch 2 taken 1502 times.
|
1544 | if ( m_fstransient->iteration == 1 ) { |
| 869 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 24 times.
|
42 | if ( r_instance.useALEformulation ) { |
| 870 | 18 | m_tau = 1./(m_fstransient->timeStep); | |
| 871 | 18 | m_cycl = 1; // For movemesh case | |
| 872 | } | ||
| 873 | |||
| 874 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 38 times.
|
42 | if (m_lumpedModelBC) { |
| 875 | // m_labelFluxListLumpedModelBC computation with LumpedModelBC labels | ||
| 876 | 4 | m_labelFluxListLumpedModelBC.clear(); | |
| 877 |
2/2✓ Branch 1 taken 8 times.
✓ Branch 2 taken 4 times.
|
12 | for(std::size_t i=0; i<m_lumpedModelBC->size(); i++) { |
| 878 | 8 | m_labelFluxListLumpedModelBC.push_back(r_instance.lumpedModelBCLabel[i]); | |
| 879 | } | ||
| 880 | } | ||
| 881 | |||
| 882 | 42 | m_bdfNS.initialize(m_pSolverNS->solution()); | |
| 883 | 42 | m_bdfNS.extrapolate(m_solExtrapolate); | |
| 884 | 42 | m_pSolverNS->initExtrapol(m_solExtrapolate); | |
| 885 | |||
| 886 | // Crank-Nicolson time marching | ||
| 887 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 42 times.
|
42 | if ( Tools::equal(r_instance.theta,0.5) ) { |
| 888 | ✗ | m_CrankNicolsonNS->initialize(m_pSolverNS->solution()); | |
| 889 | ✗ | m_CrankNicolsonNS->extrapolate(m_solExtrapolate); | |
| 890 | ✗ | m_pSolverNS->updateExtrapol(m_solExtrapolate); | |
| 891 | } | ||
| 892 | |||
| 893 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42 times.
|
42 | if ( r_instance.readPreStressedDisplFromFile ) { |
| 894 | ✗ | computePreStressedDispMeshVector(); | |
| 895 | } | ||
| 896 | } | ||
| 897 | 1544 | } | |
| 898 | |||
| 899 | /***********************************************************************************/ | ||
| 900 | /***********************************************************************************/ | ||
| 901 | |||
| 902 | ✗ | int NSModel::getNstate() const | |
| 903 | { | ||
| 904 | ✗ | return m_pSolverNS->numDof(); | |
| 905 | } | ||
| 906 | |||
| 907 | /***********************************************************************************/ | ||
| 908 | /***********************************************************************************/ | ||
| 909 | |||
| 910 | ✗ | void NSModel::getState(double* & state) | |
| 911 | { | ||
| 912 | ✗ | m_pSolverNS->getSolution(state, MpiInfo::numProc(), MpiInfo::rankProc()); | |
| 913 | } | ||
| 914 | |||
| 915 | /***********************************************************************************/ | ||
| 916 | /***********************************************************************************/ | ||
| 917 | |||
| 918 | ✗ | felInt NSModel::numDof_swig() | |
| 919 | { | ||
| 920 | ✗ | return m_pSolverNS->numDof(); | |
| 921 | } | ||
| 922 | |||
| 923 | /***********************************************************************************/ | ||
| 924 | /***********************************************************************************/ | ||
| 925 | |||
| 926 | ✗ | void NSModel::getState_swig(double* data, felInt numDof) | |
| 927 | { | ||
| 928 | double * state; | ||
| 929 | ✗ | m_pSolverNS->getSolution(state, MpiInfo::numProc(), MpiInfo::rankProc()); | |
| 930 | ✗ | for (felInt i=0 ; i<numDof ; i++) { | |
| 931 | ✗ | data[i] = state[i]; | |
| 932 | } | ||
| 933 | } | ||
| 934 | |||
| 935 | /***********************************************************************************/ | ||
| 936 | /***********************************************************************************/ | ||
| 937 | |||
| 938 | ✗ | void NSModel::setState(double* & state) | |
| 939 | { | ||
| 940 | ✗ | m_pSolverNS->setSolution(state, MpiInfo::numProc(), MpiInfo::rankProc()); | |
| 941 | } | ||
| 942 | |||
| 943 | /***********************************************************************************/ | ||
| 944 | /***********************************************************************************/ | ||
| 945 | |||
| 946 | 360 | void NSModel::computeALEterms() | |
| 947 | { | ||
| 948 | // Retrieve instance | ||
| 949 | 360 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 950 | |||
| 951 | // Assembly the RHS (flagMatrixRHS = FlagMatrixRHS::only_rhs) in the previous mesh configuration: | ||
| 952 | // \dfrac{\rho}{\dt} \int_{\Omega^{n-1}} u^{n-1} v + \int_{\Omega^{n-1}} f v | ||
| 953 | 360 | m_pSolverNS->onPreviousMesh(); | |
| 954 | 360 | m_pSolverNS->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::only_rhs); | |
| 955 | |||
| 956 | //Save velocity previous mesh | ||
| 957 |
2/2✓ Branch 0 taken 160 times.
✓ Branch 1 taken 200 times.
|
360 | if ( r_instance.nonLinearFluid ) { |
| 958 | 160 | m_solPrev.zeroEntries(); | |
| 959 | 160 | m_solPrev.copyFrom(m_pSolverNS->vector()); | |
| 960 | } | ||
| 961 | |||
| 962 | //First distinguish the elastic extension and harmonic extension. | ||
| 963 |
2/2✓ Branch 0 taken 70 times.
✓ Branch 1 taken 290 times.
|
360 | if ( r_instance.useElasticExtension ) { |
| 964 | 70 | UpdateMeshByDisplacement(); | |
| 965 | } | ||
| 966 | else{ | ||
| 967 | // Solve the extension problem and move the mesh of m_pSolverNS | ||
| 968 | // according to the mesh displacement obtained | ||
| 969 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 290 times.
|
290 | if( r_instance.updateMeshByVelocity ) { |
| 970 | // the HarmonicExtention unknown is velocity | ||
| 971 | // dispalcement is computed analitycally | ||
| 972 | ✗ | UpdateMeshByVelocity(); | |
| 973 | } else { | ||
| 974 | // the HarmonicExtention unknown is displacement | ||
| 975 | // velocity is computed analitycally | ||
| 976 | 290 | UpdateMeshByDisplacement(); | |
| 977 | } | ||
| 978 | } | ||
| 979 | |||
| 980 | 360 | m_pSolverNS->onCurrentMesh(); | |
| 981 | |||
| 982 | // Assembly the Matrix (flagMatrixRHS = FlagMatrixRHS::only_matrix) in the current mesh configuration | ||
| 983 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 360 times.
|
360 | if ( r_instance.NSStabType == 1 ) |
| 984 | ✗ | m_pSolverNS->assembleFaceOrientedStabilization(); | |
| 985 | |||
| 986 | 360 | m_pSolverNS->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::only_matrix); | |
| 987 | 360 | m_pSolverNS->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::only_rhs); | |
| 988 | 360 | } | |
| 989 | |||
| 990 | /***********************************************************************************/ | ||
| 991 | /***********************************************************************************/ | ||
| 992 | |||
| 993 | ✗ | void NSModel::UpdateMeshByVelocity() | |
| 994 | { | ||
| 995 | // Retrieve instance | ||
| 996 | ✗ | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 997 | |||
| 998 | ✗ | if ( r_instance.hasInitialCondition ) { // m_solDispl_0 = d^{0} != 0 | |
| 999 | ✗ | m_solDispl_0.copyFrom(m_pSolverHE->solution()); | |
| 1000 | //else m_solDispl_0 = d^{0} == 0 defined in ExternalVec() | ||
| 1001 | } | ||
| 1002 | |||
| 1003 | ✗ | if(m_fstransient->iteration>1 ) { | |
| 1004 | ✗ | m_solDispl.copyFrom(m_pSolverHE->solution()); | |
| 1005 | } | ||
| 1006 | |||
| 1007 | ✗ | felInt displIterMax = r_instance.timeMaxCaseFile; // todo: it should be taken automatically form displacement file.case | |
| 1008 | |||
| 1009 | ✗ | if(m_fstransient->iteration == m_cycl * displIterMax +1) { // to read the file case in cycle | |
| 1010 | ✗ | m_displTimeStep = m_fstransient->timeStep ; | |
| 1011 | ✗ | m_cycl++; | |
| 1012 | } | ||
| 1013 | ✗ | PetscPrintf(MpiInfo::petscComm(),"Cycle %d, Read data at time t = %f\n", m_cycl, m_displTimeStep); | |
| 1014 | ✗ | m_pSolverHE->readDataDisplacement(m_ios, m_displTimeStep); // defined in linearProblemHarmonicExtention | |
| 1015 | |||
| 1016 | //Solve the extension problem | ||
| 1017 | ✗ | solveHarmExt(); | |
| 1018 | |||
| 1019 | // Store velocity and apply it to the NS problem. seqSolution is the velocity | ||
| 1020 | ✗ | m_pSolverHE->gatherSolution(); | |
| 1021 | ✗ | m_velMesh.copyFrom(m_pSolverHE->sequentialSolution()); | |
| 1022 | ✗ | m_velMesh.scale(-1.); | |
| 1023 | |||
| 1024 | /*In which follow compute displacement: d^n = w^n*dt + d^{n-1} | ||
| 1025 | VecAXPBY(Vec y,PetscScalar a,PetscScalar b,Vec x); -> y=a∗x+b∗y | ||
| 1026 | solution() == w^n | ||
| 1027 | m_solDispl == d^{n-1} | ||
| 1028 | when a cardiac cycle re-start we need to initilize the displacement(t==0) | ||
| 1029 | because the displacement(0) != dispalcement(timeMaxCaseFile), | ||
| 1030 | where timeMaxCaseFile is last time step*/ | ||
| 1031 | ✗ | if(m_fstransient->iteration == m_cycl * displIterMax + 1) | |
| 1032 | ✗ | m_pSolverHE->solution().axpby(1, r_instance.timeStep, m_solDispl_0); | |
| 1033 | else | ||
| 1034 | ✗ | m_pSolverHE->solution().axpby(1, r_instance.timeStep, m_solDispl); | |
| 1035 | |||
| 1036 | ✗ | m_pSolverHE->gatherSolution(); | |
| 1037 | // seqSolution() is now the displacement because in VecAXPBY we change solution() | ||
| 1038 | |||
| 1039 | ✗ | updateMesh(); | |
| 1040 | } | ||
| 1041 | |||
| 1042 | /***********************************************************************************/ | ||
| 1043 | /***********************************************************************************/ | ||
| 1044 | |||
| 1045 | 360 | void NSModel::UpdateMeshByDisplacement() | |
| 1046 | { | ||
| 1047 | // Retrieve instance | ||
| 1048 | 360 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 1049 | |||
| 1050 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 360 times.
|
360 | if ( r_instance.readDisplFromFile ) |
| 1051 | ✗ | m_pSolverHE->readDataDisplacement(m_ios, m_displTimeStep); // defined in linearProblemHarmonicExtention | |
| 1052 | |||
| 1053 | // fsi challenge | ||
| 1054 |
2/2✓ Branch 0 taken 70 times.
✓ Branch 1 taken 290 times.
|
360 | if ( r_instance.useElasticExtension ) |
| 1055 | 70 | m_dispMesho_pp.copyFrom( m_pSolverHE->solution() ); | |
| 1056 | |||
| 1057 | |||
| 1058 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 360 times.
|
360 | if ( Tools::equal(r_instance.theta,0.5) ) { |
| 1059 | ✗ | m_velMesh.copyFrom(m_dispMesho); | |
| 1060 | ✗ | m_dispMesho.copyFrom(m_dispMesh); | |
| 1061 | } else { | ||
| 1062 | // m_dispMesh = 0 at m_fstransient->iteration = 1 | ||
| 1063 | 360 | m_velMesh.copyFrom(m_dispMesh); | |
| 1064 | } | ||
| 1065 | |||
| 1066 | //Solve the extension problem | ||
| 1067 | 360 | solveHarmExt(); | |
| 1068 | 360 | m_pSolverHE->gatherSolution(); | |
| 1069 | |||
| 1070 | 360 | updateMesh(); | |
| 1071 | |||
| 1072 | |||
| 1073 | // Compute velocity and apply it to the NS problem. seqSolution is the displacement | ||
| 1074 | //VecAXPBY(Vec y,PetscScalar a,PetscScalar b,Vec x); -> y=a∗x+b∗y | ||
| 1075 | // now m_velMesh == d^{n-1} | ||
| 1076 | // and m_dispMesh == d^{n} | ||
| 1077 | 360 | m_velMesh.axpby(-m_tau,m_tau,m_dispMesh); | |
| 1078 | // Result: m_velMesh = -w^{n} = \dfrac{1}{\dt} (-d^{n} + d^{n-1}) | ||
| 1079 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 360 times.
|
360 | if ( Tools::equal(r_instance.theta,0.5) ) |
| 1080 | ✗ | m_velMesh.scale(0.5); | |
| 1081 | // Result: m_velMesh = -w^{n+1/2} = 0.5 \dfrac{1}{\dt} (-d^{n} + d^{n-2}) | ||
| 1082 | 360 | } | |
| 1083 | |||
| 1084 | /***********************************************************************************/ | ||
| 1085 | /***********************************************************************************/ | ||
| 1086 | |||
| 1087 | 360 | void NSModel::solveHarmExt() | |
| 1088 | { | ||
| 1089 | // Retrieve instance | ||
| 1090 | 360 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 1091 | |||
| 1092 | //Assemble matrix and RHS | ||
| 1093 | 360 | m_pSolverHE->assembleMatrixRHS(MpiInfo::rankProc()); | |
| 1094 | //Specific operations before solve the system | ||
| 1095 | 360 | postAssemblingMatrixRHS(1); | |
| 1096 | //Apply essential transient boundary conditions | ||
| 1097 | 360 | m_pSolverHE->finalizeEssBCTransient(); | |
| 1098 | 360 | m_pSolverHE->applyBC(r_instance.essentialBoundaryConditionsMethod, MpiInfo::rankProc()); | |
| 1099 | //Solve linear system | ||
| 1100 | 360 | m_pSolverHE->solve(MpiInfo::rankProc(), MpiInfo::numProc()); | |
| 1101 | |||
| 1102 | // fsi challenge (update displacement with increment) | ||
| 1103 |
2/2✓ Branch 0 taken 70 times.
✓ Branch 1 taken 290 times.
|
360 | if ( r_instance.useElasticExtension ) |
| 1104 | 70 | m_pSolverHE->solution().axpy(1.,m_dispMesho_pp); | |
| 1105 | 360 | } | |
| 1106 | |||
| 1107 | /***********************************************************************************/ | ||
| 1108 | /***********************************************************************************/ | ||
| 1109 | |||
| 1110 | 360 | void NSModel::updateMesh() | |
| 1111 | { | ||
| 1112 | // Retrieve instance | ||
| 1113 | 360 | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 1114 | |||
| 1115 | 360 | const int ivar = m_pSolverHE->listVariable().getVariableIdList(displacement); | |
| 1116 | 360 | const int iunk = m_pSolverHE->listUnknown().getUnknownIdList(displacement); | |
| 1117 | 360 | const int iMesh = m_pSolverHE->listVariable()[ivar].idMesh(); | |
| 1118 | |||
| 1119 | //========= | ||
| 1120 | // ALE | ||
| 1121 | //========= | ||
| 1122 | 360 | PetscPrintf(MpiInfo::petscComm(),"Update mesh\n"); | |
| 1123 | |||
| 1124 | 360 | m_dispMesh.copyFrom(m_pSolverHE->sequentialSolution()); | |
| 1125 | // m_dispMesh is the solution obtained with solveHarmExt() | ||
| 1126 | // Extract the values of m_dispMesh and put them in the std::vector m_dispMeshVector | ||
| 1127 | // according to the application ordering of the nodes | ||
| 1128 | 360 | m_dispMesh.getValues(m_numDofExtensionProblem, m_petscToGlobal.data(), m_tmpDisp.data()); | |
| 1129 | |||
| 1130 | 360 | felInt numNodes = m_pSolverHE->supportDofUnknown(iunk).listNode().size(); | |
| 1131 | felInt idof; | ||
| 1132 | 360 | felInt numCompDisp = m_pSolverHE->dimension(); | |
| 1133 | 360 | felInt numCoorMesh = m_pSolverHE->meshLocal(iMesh)->numCoor(); | |
| 1134 | |||
| 1135 |
2/2✓ Branch 0 taken 40 times.
✓ Branch 1 taken 320 times.
|
360 | if ( r_instance.FusionDof ) { |
| 1136 |
2/2✓ Branch 0 taken 40800 times.
✓ Branch 1 taken 40 times.
|
40840 | for (felInt ip = 0; ip < numNodes; ip++) { |
| 1137 | 40800 | idof = m_pSolverHE->supportDofUnknown(iunk).listPerm()[ip]; | |
| 1138 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 40800 times.
|
40800 | if ( numCompDisp < numCoorMesh ) { |
| 1139 | ✗ | m_dispMeshVector[ip*3] = m_tmpDisp[idof*2]; | |
| 1140 | ✗ | m_dispMeshVector[ip*3+1] = m_tmpDisp[idof*2+1]; | |
| 1141 | ✗ | m_dispMeshVector[ip*3+2] = 0.; | |
| 1142 | } | ||
| 1143 |
1/2✓ Branch 0 taken 40800 times.
✗ Branch 1 not taken.
|
40800 | else if ( numCompDisp == numCoorMesh ) { |
| 1144 |
2/2✓ Branch 0 taken 122400 times.
✓ Branch 1 taken 40800 times.
|
163200 | for (felInt ic = 0; ic < numCoorMesh; ++ic) |
| 1145 | 122400 | m_dispMeshVector[ip*numCoorMesh+ic] = m_tmpDisp[idof*numCoorMesh+ic]; | |
| 1146 | } | ||
| 1147 | } | ||
| 1148 | } | ||
| 1149 | else // to be modified: does not work if numCompDisp < numCoorMesh // TODO | ||
| 1150 | 320 | m_dispMeshVector = m_tmpDisp; | |
| 1151 | |||
| 1152 | //Move the mesh in the fluid problem | ||
| 1153 |
1/2✓ Branch 0 taken 360 times.
✗ Branch 1 not taken.
|
360 | if ( r_instance.NSequationFlag != 0 ) { |
| 1154 |
1/2✓ Branch 3 taken 360 times.
✗ Branch 4 not taken.
|
360 | m_pSolverNS->meshLocal(iMesh)->moveMesh(m_dispMeshVector,1.0); |
| 1155 |
4/6✓ Branch 1 taken 115 times.
✓ Branch 2 taken 245 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 115 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 360 times.
|
360 | if ( (MpiInfo::rankProc() == 0) && r_instance.hasImmersedData ) // Attention: we move also the points of the global mesh, this can impact the fluid.geo if written at each time-step |
| 1156 | ✗ | m_pSolverNS->mesh(iMesh)->moveMesh(m_dispMeshVector,1.0); | |
| 1157 |
2/2✓ Branch 0 taken 70 times.
✓ Branch 1 taken 290 times.
|
360 | if ( r_instance.useElasticExtension ) |
| 1158 | // extension computed on the deformed mesh (takes into account the size of the deformed elements) | ||
| 1159 |
1/2✓ Branch 3 taken 70 times.
✗ Branch 4 not taken.
|
70 | m_pSolverHE->meshLocal(iMesh)->moveMesh(m_dispMeshVector,1.0); |
| 1160 | } | ||
| 1161 | 360 | } | |
| 1162 | |||
| 1163 | /***********************************************************************************/ | ||
| 1164 | /***********************************************************************************/ | ||
| 1165 | |||
| 1166 | ✗ | void NSModel::computePreStressedDispMeshVector() | |
| 1167 | { | ||
| 1168 | // Retrieve instance | ||
| 1169 | ✗ | auto& r_instance = FelisceParam::instance(instanceIndex()); | |
| 1170 | |||
| 1171 | ✗ | m_pSolverHE->HarmExtSol().getValues(m_numDofExtensionProblem, m_petscToGlobal.data(), m_tmpDisp.data()); | |
| 1172 | |||
| 1173 | ✗ | const int iunk = m_pSolverHE->listUnknown().getUnknownIdList(displacement); | |
| 1174 | |||
| 1175 | ✗ | felInt numNodes = m_pSolverHE->supportDofUnknown(iunk).listNode().size(); | |
| 1176 | ✗ | if(r_instance.FusionDof) { | |
| 1177 | |||
| 1178 | // if there are e.g. RIS model, the number of dofs are re-numerated. | ||
| 1179 | ✗ | m_preStressedDispMeshVector.resize(numNodes*3, 0.); | |
| 1180 | ✗ | for (felInt ip = 0; ip < numNodes; ip++) { | |
| 1181 | ✗ | m_preStressedDispMeshVector[ip*3] = m_tmpDisp[m_pSolverHE->supportDofUnknown(iunk).listPerm()[ip]*3]; // TODO D.C. supportDofUnknown(0) is velocity? maybe is better get the id of the unknow from listVariable | |
| 1182 | ✗ | m_preStressedDispMeshVector[ip*3+1] = m_tmpDisp[m_pSolverHE->supportDofUnknown(iunk).listPerm()[ip]*3+1]; | |
| 1183 | ✗ | m_preStressedDispMeshVector[ip*3+2] = m_tmpDisp[m_pSolverHE->supportDofUnknown(iunk).listPerm()[ip]*3+2]; | |
| 1184 | } | ||
| 1185 | } else | ||
| 1186 | ✗ | m_preStressedDispMeshVector = m_tmpDisp; | |
| 1187 | } | ||
| 1188 | |||
| 1189 | /***********************************************************************************/ | ||
| 1190 | /***********************************************************************************/ | ||
| 1191 | |||
| 1192 | 542 | PetscVector& NSModel::seqResidual() | |
| 1193 | { | ||
| 1194 | 542 | return m_seqResidual; | |
| 1195 | } | ||
| 1196 | |||
| 1197 | /***********************************************************************************/ | ||
| 1198 | /***********************************************************************************/ | ||
| 1199 | |||
| 1200 | ✗ | PetscVector& NSModel::residual() | |
| 1201 | { | ||
| 1202 | ✗ | return m_residual; | |
| 1203 | } | ||
| 1204 | |||
| 1205 | /***********************************************************************************/ | ||
| 1206 | /***********************************************************************************/ | ||
| 1207 | |||
| 1208 | 6 | int& NSModel::hasLumpedMassRN() | |
| 1209 | { | ||
| 1210 | 6 | return m_hasLumpedMassRN; | |
| 1211 | } | ||
| 1212 | |||
| 1213 | /***********************************************************************************/ | ||
| 1214 | /***********************************************************************************/ | ||
| 1215 | |||
| 1216 | 18 | PetscVector& NSModel::lumpedMassRN() | |
| 1217 | { | ||
| 1218 | 18 | return m_lumpedMassRN; | |
| 1219 | } | ||
| 1220 | |||
| 1221 | /***********************************************************************************/ | ||
| 1222 | /***********************************************************************************/ | ||
| 1223 | |||
| 1224 | 8 | Immersed & NSModel::immersed() | |
| 1225 | { | ||
| 1226 | 8 | return m_immersed; | |
| 1227 | } | ||
| 1228 | |||
| 1229 | /***********************************************************************************/ | ||
| 1230 | /***********************************************************************************/ | ||
| 1231 | |||
| 1232 | ✗ | void NSModel::writeCurrentDisplacement(const int mesh_id, const std::string name) | |
| 1233 | { | ||
| 1234 | |||
| 1235 | ✗ | if ( MpiInfo::rankProc() != 0 ) | |
| 1236 | ✗ | return; | |
| 1237 | |||
| 1238 | ✗ | const int points = m_mesh[mesh_id]->numPoints(); | |
| 1239 | ✗ | const int dim = 3*points; | |
| 1240 | ✗ | double* disp = new double[dim]; | |
| 1241 | |||
| 1242 | ✗ | for (int i = 0; i < points; ++i){ | |
| 1243 | ✗ | for (int j = 0; j < 3; ++j) | |
| 1244 | ✗ | disp[3*i+j] = m_mesh[mesh_id]->listPoints()[i].coor(j) - m_mesh[mesh_id]->listReferencePoints()[i].coor(j); | |
| 1245 | } | ||
| 1246 | ✗ | m_ios[mesh_id]->writeSolution(MpiInfo::rankProc(), m_fstransient->time, m_fstransient->iteration, 1, name , disp, points, 3, true); | |
| 1247 | } | ||
| 1248 | |||
| 1249 | /***********************************************************************************/ | ||
| 1250 | /***********************************************************************************/ | ||
| 1251 | |||
| 1252 | ✗ | void NSModel::writeCurrentNormalVector() | |
| 1253 | { | ||
| 1254 | ✗ | if ( MpiInfo::rankProc() != 0 ) | |
| 1255 | ✗ | return; | |
| 1256 | |||
| 1257 | ✗ | const felInt iUnknownLag = m_pSolverNS->iUnknownLag(); | |
| 1258 | ✗ | const felInt iLagrange = m_pSolverNS->idVarLag(); | |
| 1259 | ✗ | const felInt lagMesh = m_pSolverNS->listVariable()[iLagrange].idMesh(); | |
| 1260 | ✗ | const int nodes = m_pSolverNS->supportDofUnknown(iUnknownLag).listNode().size(); | |
| 1261 | ✗ | const int dim = 3*nodes; | |
| 1262 | ✗ | double* normal = new double[dim]; | |
| 1263 | ✗ | for (int i = 0; i < nodes; ++i){ | |
| 1264 | ✗ | for (int j = 0; j < 3; ++j) | |
| 1265 | ✗ | normal[3*i+j] = m_pSolverNS->mesh(lagMesh)->listNormal(i).coor(j); | |
| 1266 | } | ||
| 1267 | ✗ | std::string name = "normal"; | |
| 1268 | ✗ | m_ios[lagMesh]->writeSolution(MpiInfo::rankProc(), m_fstransient->time, m_fstransient->iteration, 1, name , normal, nodes, 3, true); | |
| 1269 | } | ||
| 1270 | |||
| 1271 | /***********************************************************************************/ | ||
| 1272 | /************************* CVGRAPH *************************************************/ | ||
| 1273 | |||
| 1274 | #ifdef FELISCE_WITH_CVGRAPH | ||
| 1275 | ✗ | void NSModel::cvgraphNewTimeStep() | |
| 1276 | { | ||
| 1277 | ✗ | this->prepareForward(); | |
| 1278 | } | ||
| 1279 | #endif | ||
| 1280 | } | ||
| 1281 |