| Directory: | ./ |
|---|---|
| File: | Solver/linearProblemPerfectFluidRS.cpp |
| Date: | 2024-04-14 07:32:34 |
| Exec | Total | Coverage | |
|---|---|---|---|
| Lines: | 41 | 89 | 46.1% |
| Branches: | 70 | 234 | 29.9% |
| 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 "Solver/linearProblemPerfectFluidRS.hpp" | ||
| 21 | #include "DegreeOfFreedom/boundaryCondition.hpp" | ||
| 22 | |||
| 23 | namespace felisce{ | ||
| 24 | /* \brief Function to assemble the mass matrix | ||
| 25 | * Function to be called in the assemblyLoopBoundaryGeneral | ||
| 26 | */ | ||
| 27 | ✗ | void LinearProblemPerfectFluidRS::massMatrixComputer(felInt ielSupportDof) { | |
| 28 | ✗ | this->m_elementMatBD[0]->zero(); | |
| 29 | ✗ | this->m_elementMatBD[0]->phi_i_phi_j(/*coef*/1.,*m_feIopCurv,/*iblock*/0,/*iblock*/0,/*numComponent*/1.); | |
| 30 | ✗ | this->setValueMatrixBD(ielSupportDof); | |
| 31 | } | ||
| 32 | |||
| 33 | /* \brief Function to assemble the laplacian matrix | ||
| 34 | * Function to be called in the assemblyLoopBoundaryGeneral | ||
| 35 | */ | ||
| 36 | ✗ | void LinearProblemPerfectFluidRS::laplacianMatrixComputer(felInt ielSupportDof) { | |
| 37 | ✗ | for ( felInt efComp(0); efComp < this->dimension(); ++efComp ) { | |
| 38 | ✗ | std::stringstream name; | |
| 39 | ✗ | name<<"normalVector" << efComp; | |
| 40 | ✗ | m_normalField.setValue(this->m_seqVecs.Get(name.str()), *m_feIopCurv , ielSupportDof, /*m_idVarIop*/0, this->m_ao, this->dof(), | |
| 41 | /*component where to read in inputvector*/0, /*where to save in elemfield*/efComp); | ||
| 42 | } | ||
| 43 | ✗ | m_curv.update(m_normalField, m_feIopCurv->covMetric[0], m_feIopCurv->m_covBasis[0]); | |
| 44 | ✗ | this->m_elementMatBD[0]->zero(); | |
| 45 | ✗ | this->m_elementMatBD[0]->sGrad_psi_j_tensor_sGrad_phi_i_for_scalar(/*coef*/ 1, *m_feIopCurv, m_curv.invMetricTensor(),/*idBlock*/0 ); | |
| 46 | ✗ | this->setValueMatrixBD(ielSupportDof); | |
| 47 | } | ||
| 48 | void | ||
| 49 | ✗ | LinearProblemPerfectFluidRS::initPerETLAP() { | |
| 50 | ✗ | this->initPerETMASS(); | |
| 51 | ✗ | m_normalField.initialize( DOF_FIELD, *m_feIopCurv, m_feIopCurv->numCoor() ); | |
| 52 | } | ||
| 53 | void | ||
| 54 | ✗ | LinearProblemPerfectFluidRS::initPerETMASS() { | |
| 55 | ✗ | felInt numDofTotal = 0; | |
| 56 | //pay attention this numDofTotal is bigger than expected | ||
| 57 | ✗ | for (std::size_t iFe = 0; iFe < this->m_listCurvilinearFiniteElement.size(); iFe++) {//this loop is on variables while it should be on unknown | |
| 58 | ✗ | numDofTotal += this->m_listCurvilinearFiniteElement[iFe]->numDof()*this->m_listVariable[iFe].numComponent(); | |
| 59 | } | ||
| 60 | ✗ | m_globPosColumn.resize(numDofTotal); m_globPosRow.resize(numDofTotal); m_matrixValues.resize(numDofTotal*numDofTotal); | |
| 61 | |||
| 62 | ✗ | m_feIopCurv = this->m_listCurvilinearFiniteElement[ 0 /* m_idVarIop */ ]; | |
| 63 | } | ||
| 64 | void | ||
| 65 | ✗ | LinearProblemPerfectFluidRS::updateFE(const std::vector<Point*>& elemPoint, const std::vector<int>&) { | |
| 66 | ✗ | m_feIopCurv->updateMeasNormal(0, elemPoint); | |
| 67 | } | ||
| 68 | |||
| 69 | /*! \brief DofBoundary class is built | ||
| 70 | * This function calls the methods to correctly initialize the DofBoundary class. | ||
| 71 | */ | ||
| 72 | void | ||
| 73 | 4 | LinearProblemPerfectFluidRS::initializeDofBoundaryAndBD2VolMaps() { | |
| 74 | /// this function builds all the mappings | ||
| 75 | /// to deal with petsc objects defined at the interface | ||
| 76 | 4 | m_dofBD[0/*iBD*/].initialize(this); | |
| 77 |
1/2✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
|
4 | this->dofBD(/*iBD*/0).buildListOfBoundaryPetscDofs(this, FelisceParam::instance().labelIOPMesh, this->idUnknownIop() , /*iComp*/0); |
| 78 | 4 | m_dofBD[0/*iBD*/].buildBoundaryVolumeMapping( m_iUnknownIop, 0); | |
| 79 | 4 | m_interfaceLabels=FelisceParam::instance().labelIOPMesh; | |
| 80 | 4 | } | |
| 81 | // this function is supposed to be called in the model | ||
| 82 | // and it solves the problem using the S-P operator | ||
| 83 | // we assume that the steklov is allocated | ||
| 84 | // and that velocity interface has been already read | ||
| 85 | // from NS | ||
| 86 | void | ||
| 87 | 152 | LinearProblemPerfectFluidRS::solvePFUsingSteklov(FelisceTransient::Pointer fs, int nfp) | |
| 88 | { | ||
| 89 | // initialization and computation | ||
| 90 |
2/4✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 152 times.
✗ Branch 6 not taken.
|
152 | m_seqVecs.Init("velPlusAlphaIop"); |
| 91 |
5/10✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 152 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 152 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 152 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 152 times.
✗ Branch 16 not taken.
|
152 | m_seqVecs.Get ("velPlusAlphaIop").duplicateFrom(m_seqVecs.Get("velocityInterface")); |
| 92 |
5/10✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 152 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 152 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 152 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 152 times.
✗ Branch 16 not taken.
|
152 | m_seqVecs.Get ("velPlusAlphaIop").copyFrom(m_seqVecs.Get("velocityInterface")); |
| 93 |
5/10✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 152 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 152 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 152 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 152 times.
✗ Branch 16 not taken.
|
152 | m_seqVecs.Get ("velPlusAlphaIop").axpy(m_relaxParam,m_seqVecs.Get("currentIop")); // this iop is the last available approximation of the iop |
| 94 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 152 times.
|
152 | if ( FelisceParam::instance().porousParam < 0 ) { |
| 95 | ✗ | m_seqVecs.Get ("velPlusAlphaIop").axpy(-1,m_seqVecs.Get("velocityInterfaceOld")); | |
| 96 | } | ||
| 97 | // We could have read it directly from the solution(), however this is true only if you do fixed point iterations, we have implemented also more sophisticated techniques that further update | ||
| 98 | // the current iop | ||
| 99 | |||
| 100 |
3/6✓ Branch 4 taken 152 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 152 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 152 times.
✗ Branch 11 not taken.
|
152 | this->applySteklov(m_seqVecs.Get("velPlusAlphaIop"), this->solution(),fs,nfp); |
| 101 | |||
| 102 | // destruction | ||
| 103 |
2/4✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 152 times.
✗ Branch 6 not taken.
|
152 | m_seqVecs.erase("velPlusAlphaIop"); |
| 104 | 152 | } | |
| 105 | |||
| 106 | /*! \brief Function to apply boundary conditions | ||
| 107 | * Very important function when working with the Steklov operator. | ||
| 108 | * | ||
| 109 | * Here the bounadry conditions are detailed, depending on the | ||
| 110 | * datafile, we use different couplings. | ||
| 111 | */ | ||
| 112 | void | ||
| 113 | 19680 | LinearProblemPerfectFluidRS::userElementComputeNaturalBoundaryCondition( const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/,felInt& iel, int label) { | |
| 114 | // declaration of a pointer to BC | ||
| 115 | BoundaryCondition* BC; | ||
| 116 | //! We check among the Neumann boundary conditions if the labels in FelisceParam.labelIOPMesh are here. | ||
| 117 |
2/2✓ Branch 1 taken 19680 times.
✓ Branch 2 taken 19680 times.
|
39360 | for (std::size_t iInterface=0; iInterface<m_boundaryConditionList.numNeumannBoundaryCondition(); iInterface++ ) { |
| 118 | 19680 | BC = m_boundaryConditionList.Neumann(iInterface); | |
| 119 | // if the current label is from neumann and if it is in the interface | ||
| 120 | // and if it is of type EnsightFile | ||
| 121 |
2/4✓ Branch 7 taken 19680 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 13120 times.
|
32800 | if ( std::find( BC->listLabel().begin(), BC->listLabel().end(),label) != BC->listLabel().end() && |
| 122 |
7/12✓ Branch 0 taken 13120 times.
✓ Branch 1 taken 6560 times.
✓ Branch 3 taken 13120 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 13120 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 13120 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 13120 times.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 19680 times.
|
32800 | std::find( FelisceParam::instance().labelIOPMesh.begin(), FelisceParam::instance().labelIOPMesh.end(),label) != FelisceParam::instance().labelIOPMesh.end() && |
| 123 | ✗ | BC->typeValueBC() == EnsightFile | |
| 124 | ) { | ||
| 125 | ✗ | CurvilinearFiniteElement *cfeiop = m_listCurvilinearFiniteElement[m_idVarIop]; | |
| 126 | ✗ | if ( m_useSteklovData ) { | |
| 127 | ✗ | m_elemFieldNeumann[iInterface].setValue( m_seqVecs.Get("dataSteklov"), *cfeiop , iel, m_idVarIop, m_ao, dof()); | |
| 128 | } | ||
| 129 | else { | ||
| 130 | //! We add the coupling condition (same for Neumann and Robin-coupling | ||
| 131 | ✗ | m_elemFieldNeumann[iInterface].setValue( m_seqVecs.Get("velocityInterface"), *cfeiop, iel, m_idVarIop, m_ao, dof()); | |
| 132 | |||
| 133 | //! In the case of a porous media the above line is enough | ||
| 134 | //! With a perfect fluid we need to subtract the old value: they are already divided by the timeStep: we are approximating the acceleration. | ||
| 135 | ✗ | if ( FelisceParam::instance().porousParam < 0 ) { | |
| 136 | ✗ | ElementField oldVelInterface; | |
| 137 | ✗ | oldVelInterface.initialize(DOF_FIELD,*cfeiop,1); | |
| 138 | ✗ | oldVelInterface.setValue( m_seqVecs.Get("velocityInterfaceOld"), *cfeiop, iel, m_idVarIop, m_ao, dof()); | |
| 139 | ✗ | for ( int iLocDof(0); iLocDof<cfeiop->numDof(); ++iLocDof) { | |
| 140 | ✗ | m_elemFieldNeumann[iInterface].val(0,iLocDof) = m_elemFieldNeumann[iInterface].val(0,iLocDof) - oldVelInterface.val( 0, iLocDof ); | |
| 141 | } | ||
| 142 | } | ||
| 143 | } | ||
| 144 | } | ||
| 145 | } | ||
| 146 | //! Then we check among the Robin boundary conditions | ||
| 147 |
2/2✓ Branch 1 taken 19680 times.
✓ Branch 2 taken 19680 times.
|
39360 | for (std::size_t iInterface=0; iInterface<m_boundaryConditionList.numRobinBoundaryCondition(); iInterface++ ) { |
| 148 | 19680 | BC = m_boundaryConditionList.Robin(iInterface); | |
| 149 | // if the current label is from robin and if it is in the interface | ||
| 150 | // and if it is and EnsightFile type | ||
| 151 |
2/4✓ Branch 7 taken 19680 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6560 times.
✗ Branch 11 not taken.
|
26240 | if ( std::find( BC->listLabel().begin(), BC->listLabel().end(),label) != BC->listLabel().end() && |
| 152 |
8/12✓ Branch 0 taken 6560 times.
✓ Branch 1 taken 13120 times.
✓ Branch 3 taken 6560 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 6560 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 6560 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 6560 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 6560 times.
✓ Branch 19 taken 13120 times.
|
26240 | std::find( FelisceParam::instance().labelIOPMesh.begin(), FelisceParam::instance().labelIOPMesh.end(),label) != FelisceParam::instance().labelIOPMesh.end() && |
| 153 |
1/2✓ Branch 1 taken 6560 times.
✗ Branch 2 not taken.
|
6560 | BC->typeValueBC() == EnsightFile |
| 154 | ) { | ||
| 155 | //! the FP.alphaRobin is overwritten with the relaxation parameter. | ||
| 156 | 6560 | FelisceParam::instance().alphaRobin[iInterface] = m_relaxParam; | |
| 157 | |||
| 158 | 6560 | CurvilinearFiniteElement *cfeiop = m_listCurvilinearFiniteElement[m_idVarIop]; | |
| 159 |
1/2✓ Branch 0 taken 6560 times.
✗ Branch 1 not taken.
|
6560 | if ( m_useSteklovData ) { |
| 160 |
3/6✓ Branch 4 taken 6560 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6560 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6560 times.
✗ Branch 11 not taken.
|
6560 | m_elemFieldRobin[iInterface].setValue( m_seqVecs.Get("dataSteklov") , *cfeiop , iel, m_idVarIop, m_ao, dof()); |
| 161 | } | ||
| 162 | else { | ||
| 163 | //! We add the coupling condition (same for Neumann and Robin-coupling | ||
| 164 | ✗ | m_elemFieldRobin[iInterface].setValue( m_seqVecs.Get("velocityInterface") , *cfeiop , iel, m_idVarIop, m_ao, dof()); | |
| 165 | |||
| 166 | //! In the case of a porous media the above line is enough | ||
| 167 | //! With a perfect fluid we need to subtract the old value: they are already divided by the timeStep: we are approximating the acceleration. | ||
| 168 | ✗ | if ( FelisceParam::instance().porousParam < 0 ) { | |
| 169 | ✗ | ElementField oldVelInterface; | |
| 170 | ✗ | oldVelInterface.initialize(DOF_FIELD,*cfeiop,1); | |
| 171 | ✗ | oldVelInterface.setValue( m_seqVecs.Get("velocityInterfaceOld"), *cfeiop, iel, m_idVarIop, m_ao, dof()); | |
| 172 | ✗ | for ( int iLocDof(0); iLocDof<cfeiop->numDof(); ++iLocDof) { | |
| 173 | ✗ | m_elemFieldRobin[iInterface].val(0,iLocDof) += - oldVelInterface.val( 0, iLocDof ); | |
| 174 | } | ||
| 175 | } | ||
| 176 | |||
| 177 | //! This part differ from the Neumann coupling: the rhs contribution due to the relaxation is added | ||
| 178 | ✗ | ElementField estimatedIOP; | |
| 179 | ✗ | estimatedIOP.initialize(DOF_FIELD,*cfeiop,1); | |
| 180 | ✗ | estimatedIOP.setValue(m_seqVecs.Get("currentIop"),*cfeiop,iel,m_idVarIop,m_ao,dof()); | |
| 181 | ✗ | for ( int iLocDof(0); iLocDof<cfeiop->numDof(); ++iLocDof) { | |
| 182 | ✗ | m_elemFieldRobin[iInterface].val(0,iLocDof) += m_relaxParam * estimatedIOP.val( 0, iLocDof ); | |
| 183 | } | ||
| 184 | } | ||
| 185 | } | ||
| 186 | } | ||
| 187 | 19680 | } | |
| 188 | void | ||
| 189 | 4 | LinearProblemPerfectFluidRS::exportNormalField() { | |
| 190 |
1/2✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
4 | PetscVector tmpVec = this->dofBD(/*iBD*/0).allocateBoundaryVector(DofBoundary::parallel); |
| 191 |
2/2✓ Branch 1 taken 8 times.
✓ Branch 2 taken 4 times.
|
12 | for ( felInt iComp(0); iComp < this->dimension(); ++iComp ) { |
| 192 |
2/4✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
|
8 | std::stringstream filenameBin,namevector; |
| 193 |
2/4✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
|
8 | filenameBin<< "normal_"<<iComp; |
| 194 |
2/4✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
|
8 | namevector << "normalVector" << iComp; |
| 195 |
3/6✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
|
8 | this->dofBD(/*iBD*/0).restrictOnBoundary(this->m_seqVecs.Get(namevector.str()), tmpVec); |
| 196 |
1/2✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
8 | if ( this->dofBD(/*iBD*/0).hasDofsOnBoundary() ) |
| 197 |
5/10✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 17 taken 8 times.
✗ Branch 18 not taken.
|
8 | tmpVec.saveInBinaryFormat( this->m_dofBD[0/*iBD*/].comm(), filenameBin.str().c_str(), FelisceParam::instance().resultDir ); |
| 198 | 8 | } | |
| 199 | 4 | } | |
| 200 | } | ||
| 201 |