| Line |
Branch |
Exec |
Source |
| 1 |
|
|
// ______ ______ _ _ _____ ______ |
| 2 |
|
|
// | ____| ____| | (_)/ ____| | ____| |
| 3 |
|
|
// | |__ | |__ | | _| (___ ___| |__ |
| 4 |
|
|
// | __| | __| | | | |\___ \ / __| __| |
| 5 |
|
|
// | | | |____| |____| |____) | (__| |____ |
| 6 |
|
|
// |_| |______|______|_|_____/ \___|______| |
| 7 |
|
|
// Finite Elements for Life Sciences and Engineering |
| 8 |
|
|
// |
| 9 |
|
|
// License: LGL2.1 License |
| 10 |
|
|
// FELiScE default license: LICENSE in root folder |
| 11 |
|
|
// |
| 12 |
|
|
// Main authors: |
| 13 |
|
|
// |
| 14 |
|
|
|
| 15 |
|
|
// System includes |
| 16 |
|
|
|
| 17 |
|
|
// External includes |
| 18 |
|
|
|
| 19 |
|
|
// Project includes |
| 20 |
|
|
#include "Model/stokesLinearElasticityCoupledModel.hpp" |
| 21 |
|
|
#include "InputOutput/io.hpp" |
| 22 |
|
|
|
| 23 |
|
|
namespace felisce { |
| 24 |
|
|
// ====================================================================== |
| 25 |
|
|
// the public methods |
| 26 |
|
|
// ====================================================================== |
| 27 |
|
|
|
| 28 |
|
|
// constructor: it first call constructor of mother class |
| 29 |
|
✗ |
StokesLinearElasticityCoupledModel::StokesLinearElasticityCoupledModel():LinearElasticityModel(),m_timeChecker(nullptr) { |
| 30 |
|
|
// change the name of the model |
| 31 |
|
✗ |
m_name="Linear Elasticity coupled with Stokes equations"; |
| 32 |
|
|
// fixing the id of the stokes problem |
| 33 |
|
|
// this should be done automatically based on |
| 34 |
|
|
// the names of the linearProblems or on |
| 35 |
|
|
// same check on the type. |
| 36 |
|
✗ |
m_idLE=0; |
| 37 |
|
✗ |
m_idStokes=1; |
| 38 |
|
✗ |
m_interfaceLabels=FelisceParam::instance().fsiInterfaceLabel; |
| 39 |
|
|
} |
| 40 |
|
|
|
| 41 |
|
|
|
| 42 |
|
|
// this function override LinearElasticityModel::initializeDerivedModel() |
| 43 |
|
|
// it is then call in initializeModel in the Model class |
| 44 |
|
|
void |
| 45 |
|
✗ |
StokesLinearElasticityCoupledModel::initializeDerivedModel() { |
| 46 |
|
|
// first a call to the ovverriden method to initialize the mother class |
| 47 |
|
✗ |
LinearElasticityModel::initializeDerivedModel(); |
| 48 |
|
|
|
| 49 |
|
|
// Two static casts to initialize the pointers |
| 50 |
|
|
// they are usufull to call functions defined in derived linear problems, |
| 51 |
|
|
// but not declared in linearProblem class. |
| 52 |
|
✗ |
m_lpbLE = static_cast<LinearProblemLinearElasticity*> ( m_linearProblem[ m_idLE ]); |
| 53 |
|
✗ |
m_lpbStokes = static_cast<LinearProblemNSRS*> ( m_linearProblem[ m_idStokes ] ); |
| 54 |
|
|
|
| 55 |
|
|
// Setting a flag in LE linProb to true to switch several functionalities on |
| 56 |
|
✗ |
m_lpbLE->coupled(true); |
| 57 |
|
✗ |
m_lpbLE->initPetscVecsForCoupling(); |
| 58 |
|
|
// The dofBoundary object of each linear problem is initialized. |
| 59 |
|
✗ |
m_lpbLE ->initializeDofBoundaryAndBD2VolMaps(); |
| 60 |
|
✗ |
m_lpbStokes ->initializeDofBoundaryAndBD2VolMaps(); |
| 61 |
|
|
// We build the std::unordered_map the store the correspondence labels of the two meshes |
| 62 |
|
|
// it needs to be called after the buildListOfBoundaryPetscDofs |
| 63 |
|
✗ |
this->defineMapLabelInterface(); |
| 64 |
|
|
|
| 65 |
|
|
// interface labels are printed on the video to check everything is all right |
| 66 |
|
✗ |
if ( FelisceParam::verbose() > 1 ) { |
| 67 |
|
✗ |
std::stringstream labelsStokes, labelsLE; |
| 68 |
|
✗ |
labelsStokes <<"Interface labels of the Stokes mesh : "; |
| 69 |
|
✗ |
labelsLE <<"Interface labels of the Linear Elasticity mesh : "; |
| 70 |
|
✗ |
for(auto it(m_labStokes2LE.begin()); it != m_labStokes2LE.end(); ++it ) { |
| 71 |
|
✗ |
labelsStokes << it->first << '\t'; |
| 72 |
|
✗ |
labelsLE << it->second << '\t'; |
| 73 |
|
|
} |
| 74 |
|
✗ |
labelsStokes << std::endl; |
| 75 |
|
✗ |
labelsLE << std::endl; |
| 76 |
|
|
|
| 77 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"%s", labelsLE.str().c_str()); |
| 78 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"%s", labelsStokes.str().c_str()); |
| 79 |
|
|
} |
| 80 |
|
|
} |
| 81 |
|
|
|
| 82 |
|
|
// Function to initialize the std::unordered_map m_labStokes2LE |
| 83 |
|
|
// In this model we make the assumption that this std::unordered_map is the identity, however be build it because of compatibility reason. |
| 84 |
|
|
// |
| 85 |
|
|
// The side effect of this function, and the real reason why so many lines are needed for this, is to check that the different mapping were correctly built in the previous function |
| 86 |
|
|
// the name of the function could be something like initMapLabelInterfaceAndCheckEverything |
| 87 |
|
|
void |
| 88 |
|
✗ |
StokesLinearElasticityCoupledModel::defineMapLabelInterface() { |
| 89 |
|
|
//saving the std::unordered_map that associates the first point (in the lexicographical order) |
| 90 |
|
|
//to each label |
| 91 |
|
✗ |
std::map<int, Point > mapStokes = m_lpbStokes->dofBD(/*iBD*/0).getMapLab2FirstPoint(); |
| 92 |
|
✗ |
std::map<int, Point > mapLE = m_lpbLE->dofBD(/*iBD*/0).getMapLab2FirstPoint(); |
| 93 |
|
|
|
| 94 |
|
|
// Checking that the number of labels from the surfaces are the same |
| 95 |
|
✗ |
std::stringstream errLine; |
| 96 |
|
✗ |
errLine<<" Different number of boundary labels describing the interfaces: "<<mapStokes.size()<< "!=" <<mapLE.size()<<std::endl; |
| 97 |
|
✗ |
FEL_CHECK(mapStokes.size()==mapLE.size(),errLine.str().c_str()); |
| 98 |
|
|
|
| 99 |
|
|
#ifndef NDEBUG |
| 100 |
|
✗ |
int counterOfNotPairedPoints(0); |
| 101 |
|
|
|
| 102 |
|
|
// For readibility, really waiting for c++11 |
| 103 |
|
|
typedef std::map<int, std::vector< Point > > map_type; |
| 104 |
|
|
|
| 105 |
|
✗ |
map_type allPointsStokes = m_lpbStokes->dofBD(/*iBD*/0).lab2AllPoints(); |
| 106 |
|
✗ |
map_type allPointsLE = m_lpbLE->dofBD(/*iBD*/0).lab2AllPoints(); |
| 107 |
|
✗ |
auto listOfPointPerLabelIteratorStokes = allPointsStokes.begin(); |
| 108 |
|
✗ |
auto listOfPointPerLabelIteratorLE = allPointsLE.begin(); |
| 109 |
|
|
|
| 110 |
|
✗ |
for ( ; /* no initialization needed */ |
| 111 |
|
✗ |
listOfPointPerLabelIteratorStokes != allPointsStokes.end() && listOfPointPerLabelIteratorLE != allPointsLE.end(); |
| 112 |
|
✗ |
++listOfPointPerLabelIteratorStokes, ++listOfPointPerLabelIteratorLE ) { |
| 113 |
|
|
|
| 114 |
|
✗ |
std::size_t numPointPerLabelStokes = listOfPointPerLabelIteratorStokes->second.size(); |
| 115 |
|
✗ |
std::size_t numPointPerLabelLE = listOfPointPerLabelIteratorLE->second.size(); |
| 116 |
|
|
|
| 117 |
|
✗ |
FEL_CHECK( numPointPerLabelStokes == numPointPerLabelLE, "Different number of points at the interface for current label\n"); |
| 118 |
|
|
|
| 119 |
|
✗ |
auto pointStokesIterator = listOfPointPerLabelIteratorStokes->second.begin(); |
| 120 |
|
✗ |
auto pointLEIterator = listOfPointPerLabelIteratorLE->second.begin(); |
| 121 |
|
|
|
| 122 |
|
✗ |
for( ; /* no initialization needed */ |
| 123 |
|
✗ |
pointStokesIterator != listOfPointPerLabelIteratorStokes->second.end() && pointLEIterator != listOfPointPerLabelIteratorLE->second.end(); |
| 124 |
|
✗ |
++pointStokesIterator, ++pointLEIterator) { |
| 125 |
|
|
|
| 126 |
|
✗ |
if ( *pointStokesIterator != *pointLEIterator ) { |
| 127 |
|
✗ |
++counterOfNotPairedPoints; |
| 128 |
|
|
} |
| 129 |
|
|
} |
| 130 |
|
✗ |
if ( counterOfNotPairedPoints > 0 ) { |
| 131 |
|
✗ |
std::stringstream aus; |
| 132 |
|
✗ |
aus<<"checking the interface: "<<counterOfNotPairedPoints<<" points differ, out of <<"<<listOfPointPerLabelIteratorLE->second.size()<<std::endl; |
| 133 |
|
✗ |
FEL_ERROR(aus.str().c_str()); |
| 134 |
|
|
} |
| 135 |
|
✗ |
MPI_Barrier(MpiInfo::petscComm()); //We wait for all the procs |
| 136 |
|
|
|
| 137 |
|
✗ |
std::stringstream msg; |
| 138 |
|
✗ |
msg<<"All right in the mapping: there is a one to one correspondance between all the nodes at the interface"<<std::endl; |
| 139 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(), "%s",msg.str().c_str()); |
| 140 |
|
|
} |
| 141 |
|
|
#endif |
| 142 |
|
|
|
| 143 |
|
|
// for each label of the interface in the Stokes mesh we check the first point |
| 144 |
|
✗ |
for(auto firstPointPerLabelStokesIterator = mapStokes.begin(); firstPointPerLabelStokesIterator != mapStokes.end(); ++firstPointPerLabelStokesIterator) { |
| 145 |
|
|
// we look for corresponding the label in the other mesh |
| 146 |
|
✗ |
bool Found=false; |
| 147 |
|
|
// for each label of the LE mesh |
| 148 |
|
✗ |
for(auto firstPointPerLabelLEIterator = mapLE.begin(); firstPointPerLabelLEIterator != mapLE.end() && !Found; ++firstPointPerLabelLEIterator) { |
| 149 |
|
|
// if the first point is the same for this couple of labels |
| 150 |
|
✗ |
if ( firstPointPerLabelStokesIterator->second == firstPointPerLabelLEIterator->second ) { |
| 151 |
|
|
// then we found the corresponding label |
| 152 |
|
✗ |
Found=true; |
| 153 |
|
|
// we fill the std::unordered_map |
| 154 |
|
✗ |
m_labStokes2LE [ firstPointPerLabelStokesIterator->first ] = firstPointPerLabelLEIterator->first; // <-- key line of the function |
| 155 |
|
|
} |
| 156 |
|
|
} |
| 157 |
|
|
// if we arrived here and we did not find a label it means that ... |
| 158 |
|
✗ |
FEL_CHECK(Found," one label of the interface of the Stokes mesh is not present in the LE mesh"); |
| 159 |
|
|
} |
| 160 |
|
|
} |
| 161 |
|
|
|
| 162 |
|
|
|
| 163 |
|
|
// Forward function: |
| 164 |
|
|
// it overrides NSSimplifiedFSI function, this time we really want to override |
| 165 |
|
|
// and we are not going to call the forward method of the mother class inside |
| 166 |
|
|
void |
| 167 |
|
✗ |
StokesLinearElasticityCoupledModel::forward() { |
| 168 |
|
|
|
| 169 |
|
|
//! First several objects and vectors that will be used through |
| 170 |
|
|
//! the different iterations are initialized/computed |
| 171 |
|
|
|
| 172 |
|
✗ |
if ( m_fstransient->iteration == 0 ) { |
| 173 |
|
|
|
| 174 |
|
|
//(1) |
| 175 |
|
|
//! Chronometers initializations |
| 176 |
|
✗ |
m_timeChecker = felisce::make_shared<ChronoInstance>(); |
| 177 |
|
✗ |
m_t1 = felisce::make_shared<ChronoInstance>(); |
| 178 |
|
✗ |
m_t2 = felisce::make_shared<ChronoInstance>(); |
| 179 |
|
✗ |
m_lpbStokes->initChronoRS(); |
| 180 |
|
|
//Total time of forward including all kind of overheads |
| 181 |
|
✗ |
m_timeChecker->start(); |
| 182 |
|
|
|
| 183 |
|
|
//(2) |
| 184 |
|
|
//! The vectors stores in m_vecs and in m_petscVecs are initialized! |
| 185 |
|
✗ |
m_t1->start(); |
| 186 |
|
✗ |
m_lpbStokes->initPetscVectors(); |
| 187 |
|
✗ |
m_t1->stop(); |
| 188 |
|
✗ |
m_t2->start(); |
| 189 |
|
✗ |
m_lpbLE->initPetscVectors(); |
| 190 |
|
|
|
| 191 |
|
|
//(3) |
| 192 |
|
|
//! The normal field is computed in m_lpbLE and then sent to the stokes. |
| 193 |
|
✗ |
m_lpbLE->computeNormalField(m_interfaceLabels,m_lpbLE->iDisplacement()); |
| 194 |
|
✗ |
m_t2->stop(); |
| 195 |
|
|
//! In the first time & fixed point iteration the P1 normal-field is exported from LE and read into Stokes |
| 196 |
|
|
//! differently from the linearProblemPerfectFluid where we do not have a std::vector variable to save |
| 197 |
|
|
//! the normal, here we have all three components of the normalField into another normalField. |
| 198 |
|
✗ |
m_t1->start(); |
| 199 |
|
✗ |
m_lpbStokes->readerInterface(m_lpbLE->dofBD(/*iBD*/0).listOfBoundaryPetscDofs(), |
| 200 |
|
✗ |
m_lpbLE->get(LinearProblem::sequential,"normalField"), |
| 201 |
|
✗ |
m_labStokes2LE, |
| 202 |
|
|
"normalField", |
| 203 |
|
✗ |
m_lpbStokes->iUnknownVel(), |
| 204 |
|
|
0/*iUnkComp0, where to read*/ |
| 205 |
|
|
); |
| 206 |
|
|
//(4) |
| 207 |
|
|
// Steklov Part! |
| 208 |
|
✗ |
if ( FelisceParam::instance().computeSteklov ) { |
| 209 |
|
✗ |
if( FelisceParam::instance().useRoSteklov ) { |
| 210 |
|
✗ |
const int ivar = m_lpbStokes->listVariable().getVariableIdList(velocity); |
| 211 |
|
✗ |
const int imesh = m_lpbStokes->listVariable()[ivar].idMesh(); |
| 212 |
|
✗ |
m_lpbStokes->getRings(); |
| 213 |
|
✗ |
m_lpbStokes->assembleLowRankSteklov(imesh); |
| 214 |
|
✗ |
m_lpbStokes->exportAllEig(imesh); |
| 215 |
|
|
} else { |
| 216 |
|
✗ |
m_lpbStokes->assembleFullRankSteklov(); |
| 217 |
|
|
} |
| 218 |
|
|
// we clean the RHS: we are going to create a new one in solveStokes |
| 219 |
|
✗ |
m_lpbStokes->clearMatrixRHS(FlagMatrixRHS::only_rhs); |
| 220 |
|
|
} |
| 221 |
|
✗ |
m_t1->stop(); |
| 222 |
|
|
|
| 223 |
|
✗ |
m_timeChecker->stop(); |
| 224 |
|
|
} |
| 225 |
|
|
|
| 226 |
|
✗ |
m_timeChecker->start(); |
| 227 |
|
|
|
| 228 |
|
|
//we do not count the time here since it is only post-processing |
| 229 |
|
|
//std::cout<<"PrepareForwardStokes: "<<m_t1->diff()<<std::endl; |
| 230 |
|
✗ |
if ( m_fstransient->iteration == 0 ) { |
| 231 |
|
|
|
| 232 |
|
✗ |
m_t1->start(); |
| 233 |
|
✗ |
this->prepareForwardStokes(); // It is important to prepare before the stokes and then LE because in stokes we check for iteration == 0 |
| 234 |
|
✗ |
m_t1->stop(); |
| 235 |
|
|
|
| 236 |
|
✗ |
m_t2->start(); |
| 237 |
|
✗ |
this->prepareForward(FlagMatrixRHS::matrix_and_rhs); // We use the prepare forward of the LEModel, we have updateTime here! There is a clean matrixrhs here! |
| 238 |
|
✗ |
m_lpbLE->assembleMatrixRHS(MpiInfo::rankProc(),FlagMatrixRHS::only_matrix); |
| 239 |
|
✗ |
m_lpbLE->addMatrixRHS(); |
| 240 |
|
✗ |
m_t2->stop(); |
| 241 |
|
|
} else { |
| 242 |
|
✗ |
m_t2->start(); |
| 243 |
|
✗ |
this->prepareForward(FlagMatrixRHS::only_rhs); // We use the prepare forward of the LEModel, we have updateTime here! There is a clean matrixrhs here! |
| 244 |
|
✗ |
m_t2->stop(); |
| 245 |
|
|
} |
| 246 |
|
|
|
| 247 |
|
|
// initialize the test quantity to one. |
| 248 |
|
|
// and the corresponding tollerance from the data file |
| 249 |
|
✗ |
double testQuantity(1.0), toll(FelisceParam::instance().domainDecompositionToll); |
| 250 |
|
|
|
| 251 |
|
|
// this two integers store the current number of sub-iterations |
| 252 |
|
|
// and of gmres iterations for the current time step |
| 253 |
|
✗ |
felInt cIteration(0), cGmresItLE(0), cGmresItStokes(0); |
| 254 |
|
|
|
| 255 |
|
|
// subiteration of domain decomposition |
| 256 |
|
✗ |
while ( testQuantity > toll && cIteration < FelisceParam::instance().domainDecompositionMaxIt) { |
| 257 |
|
|
// update the it counter and display a nice banner |
| 258 |
|
✗ |
cIteration++; |
| 259 |
|
✗ |
this->displayBannerIterationInit(cIteration); |
| 260 |
|
|
|
| 261 |
|
|
/*! LinearElasticity -- Solution Phase |
| 262 |
|
|
* |
| 263 |
|
|
* The linearElasticity problem is solved |
| 264 |
|
|
*/ |
| 265 |
|
✗ |
m_t2->start(); |
| 266 |
|
✗ |
if ( m_fstransient->iteration == 1 ) { |
| 267 |
|
✗ |
m_lpbLE->assembleMatrixRHS(MpiInfo::rankProc(),FlagMatrixRHS::only_rhs); |
| 268 |
|
✗ |
this->solveLinearElasticity(FlagMatrixRHS::matrix_and_rhs); |
| 269 |
|
|
} |
| 270 |
|
|
else { |
| 271 |
|
✗ |
m_lpbLE->assembleMatrixRHS(MpiInfo::rankProc(),FlagMatrixRHS::only_rhs); |
| 272 |
|
✗ |
this->solveLinearElasticity(FlagMatrixRHS::only_rhs); |
| 273 |
|
|
} |
| 274 |
|
✗ |
cGmresItLE += m_lpbLE->kspInterface().getNumOfIteration(); |
| 275 |
|
|
|
| 276 |
|
|
/*! LinearElasticity -- Post-processing |
| 277 |
|
|
* |
| 278 |
|
|
* Now the linear elasticity problem has been solved. |
| 279 |
|
|
* |
| 280 |
|
|
* We have to compute the velocity of the solid at the interface |
| 281 |
|
|
* The current acceleration is also computed in order to compute the velocity. |
| 282 |
|
|
*/ |
| 283 |
|
✗ |
m_lpbLE->gatherSolution(); |
| 284 |
|
✗ |
m_lpbLE->updateCurrentVelocity(); |
| 285 |
|
✗ |
m_t2->stop(); |
| 286 |
|
|
/*! Communication phase |
| 287 |
|
|
* |
| 288 |
|
|
* Now that the quantity has been computed the Stokes problem has to read it. |
| 289 |
|
|
* It is a std::vector quantity |
| 290 |
|
|
*/ |
| 291 |
|
✗ |
m_t1->start(); |
| 292 |
|
✗ |
m_lpbStokes->readerInterface( m_lpbLE->dofBD(/*iBD*/0).listOfBoundaryPetscDofs(), |
| 293 |
|
✗ |
m_lpbLE->get(LinearProblem::sequential,"velocityCurrent"), |
| 294 |
|
✗ |
m_labStokes2LE, |
| 295 |
|
|
"externalVelocity", |
| 296 |
|
✗ |
m_lpbStokes->iUnknownVel(), |
| 297 |
|
|
0/*iUnkComp0, where to read it starts from this and then it loops over the space dimension*/ |
| 298 |
|
|
); |
| 299 |
|
|
/*! Stokes -- Solution phase |
| 300 |
|
|
* |
| 301 |
|
|
* After that the data have been read we proceed and solve the Stokes problem |
| 302 |
|
|
* depending on the data file this will be done either with or without |
| 303 |
|
|
* using the Steklov operator in its full or reduced version |
| 304 |
|
|
*/ |
| 305 |
|
✗ |
this->solveStokes(cIteration); |
| 306 |
|
✗ |
cGmresItStokes += FelisceParam::instance().computeSteklov?0:m_lpbStokes->kspInterface().getNumOfIteration(); |
| 307 |
|
|
|
| 308 |
|
|
/*! Stokes -- Post-processing |
| 309 |
|
|
* |
| 310 |
|
|
* After the problem solution we need to compute the normal stress at the interface. |
| 311 |
|
|
* |
| 312 |
|
|
*/ |
| 313 |
|
✗ |
m_lpbStokes->gatherSolution(); |
| 314 |
|
✗ |
m_lpbStokes->computeResidualRS(); |
| 315 |
|
✗ |
m_t1->stop(); |
| 316 |
|
|
/*! Communication phase |
| 317 |
|
|
* |
| 318 |
|
|
* Now that the quantity has been computed the Stokes problem has to read it. |
| 319 |
|
|
* It is a std::vector quantity |
| 320 |
|
|
*/ |
| 321 |
|
✗ |
m_t2->start(); |
| 322 |
|
✗ |
m_lpbLE->accelerationPreStep(m_lpbLE->get(LinearProblem::sequential,"externalStress")); |
| 323 |
|
✗ |
m_lpbLE->readStressInterface(m_lpbStokes->dofBD(/*iBD*/0).listOfBoundaryPetscDofs(), |
| 324 |
|
✗ |
m_lpbStokes->seqResidual(), |
| 325 |
|
|
/*iUnknComp0*/0); |
| 326 |
|
|
|
| 327 |
|
|
/*! Convergence test phase |
| 328 |
|
|
* |
| 329 |
|
|
* The test quantities are the normalized L2-norms of the increments |
| 330 |
|
|
* of the two quantities exchanged at the interface: |
| 331 |
|
|
* - the velocity ( from LE to Stokes ) |
| 332 |
|
|
* - the normal stresses ( from Stokes to LE ) |
| 333 |
|
|
*/ |
| 334 |
|
✗ |
std::pair<double, double> testQuantities = m_lpbLE->computeTestQuantities(); |
| 335 |
|
✗ |
double testVelocity = testQuantities.first; |
| 336 |
|
✗ |
double testStresses = testQuantities.second; |
| 337 |
|
✗ |
testQuantity = testVelocity + testStresses; |
| 338 |
|
|
|
| 339 |
|
✗ |
m_lpbLE->accelerationPostStep( m_lpbLE->get(LinearProblem::sequential,"externalStress"), cIteration ); |
| 340 |
|
|
|
| 341 |
|
✗ |
m_t2->stop(); |
| 342 |
|
|
// we display a nice banner |
| 343 |
|
✗ |
this->displayBannerIterationEnd(cIteration,testQuantity,testVelocity,testStresses); |
| 344 |
|
|
// we clear only the RHS |
| 345 |
|
|
// it is important to clear the rhs here, because otherwise you will sum into the previous one |
| 346 |
|
|
// the former contribution |
| 347 |
|
|
// We do not clear the matrix since we want to re-use it |
| 348 |
|
✗ |
this->clearMatrixRHSOfPbs(FlagMatrixRHS::only_rhs); |
| 349 |
|
|
} |
| 350 |
|
✗ |
if ( cIteration == FelisceParam::instance().domainDecompositionMaxIt ) { |
| 351 |
|
✗ |
std::stringstream war; |
| 352 |
|
✗ |
war<<"************* Maximum number of iterations reached in "<<__FILE__<<":"<<__LINE__<<" *********************"<<std::endl; |
| 353 |
|
✗ |
FEL_WARNING(war.str().c_str()); |
| 354 |
|
|
} |
| 355 |
|
|
|
| 356 |
|
✗ |
m_timeChecker->stop(); |
| 357 |
|
|
|
| 358 |
|
|
//! we finalize the forward |
| 359 |
|
✗ |
this->finalizeForward(cIteration,cGmresItLE,cGmresItStokes); |
| 360 |
|
|
} |
| 361 |
|
|
|
| 362 |
|
✗ |
void StokesLinearElasticityCoupledModel::prepareForwardStokes() { |
| 363 |
|
✗ |
m_lpbStokes->assembleMatrixRHS(MpiInfo::rankProc()); |
| 364 |
|
|
// This is necessary since the stokes matrix is in the matrix(1) of linearProblemNS. |
| 365 |
|
|
// We have to add it all the time to the matrix(0) |
| 366 |
|
✗ |
m_lpbStokes->addMatrixRHS(); |
| 367 |
|
✗ |
m_lpbStokes->createAndCopyMatrixRHSWithoutBC(); |
| 368 |
|
|
} |
| 369 |
|
|
|
| 370 |
|
|
void |
| 371 |
|
✗ |
StokesLinearElasticityCoupledModel::writeTimeInfo() { |
| 372 |
|
|
|
| 373 |
|
|
// we convert the relaxation parameter into a std::string |
| 374 |
|
✗ |
std::stringstream aus; |
| 375 |
|
✗ |
aus<<MpiInfo::rankProc(); |
| 376 |
|
|
|
| 377 |
|
|
// we open a file in the folder before the result dir |
| 378 |
|
|
// overwriting the previous content |
| 379 |
|
|
// the information are described in the header prined here |
| 380 |
|
✗ |
std::ofstream iterationFile( std::string(FelisceParam::instance().resultDir + "timeInfo"+aus.str()+".dat").c_str() , std::ios::trunc ); |
| 381 |
|
|
iterationFile <<"TimeStep"<<'\t' |
| 382 |
|
|
<<"fixedPointIterations"<<'\t' |
| 383 |
|
|
<<"totalGmresIterationsLE"<<'\t' |
| 384 |
|
|
<<"totalGmresIterationsStokes"<<'\t' |
| 385 |
|
|
<<"gmresIterationsPerFixedPointIt"<<'\t' |
| 386 |
|
|
<<"cumulatedTime"<<'\t' |
| 387 |
|
|
<<"onlyForStokes"<<'\t' |
| 388 |
|
✗ |
<<"onlyForLE" |
| 389 |
|
✗ |
<<std::endl; |
| 390 |
|
✗ |
for ( std::size_t ts(0); ts<m_numFixedPointItPerTimeStep.size(); ++ts) { |
| 391 |
|
✗ |
iterationFile<<ts+1<<'\t' |
| 392 |
|
✗ |
<<m_numFixedPointItPerTimeStep[ts]<<'\t' |
| 393 |
|
✗ |
<<m_totNumOfGMRESItPerTimeStepLE[ts]<<'\t' |
| 394 |
|
✗ |
<<m_totNumOfGMRESItPerTimeStepStokes[ts]<<'\t' |
| 395 |
|
✗ |
<<double(m_totNumOfGMRESItPerTimeStepLE[ts] + m_totNumOfGMRESItPerTimeStepStokes[ts])/double(m_numFixedPointItPerTimeStep[ts])<<'\t' |
| 396 |
|
✗ |
<<m_cumulatedForwardTime[ts]<<'\t' |
| 397 |
|
✗ |
<<m_cumulatedStokesTime[ts]<<'\t' |
| 398 |
|
✗ |
<<m_cumulatedLETime[ts] |
| 399 |
|
✗ |
<<std::endl; |
| 400 |
|
|
} |
| 401 |
|
|
//we close the file |
| 402 |
|
✗ |
iterationFile.close(); |
| 403 |
|
|
} |
| 404 |
|
|
|
| 405 |
|
|
// ====================================================================== |
| 406 |
|
|
// now the privates methods |
| 407 |
|
|
// ====================================================================== |
| 408 |
|
|
|
| 409 |
|
|
// this function overrides |
| 410 |
|
|
// LinearElasticityModel::finalizeForward(); |
| 411 |
|
|
void |
| 412 |
|
✗ |
StokesLinearElasticityCoupledModel::finalizeForward(int cIteration, int cGmresItLE, int cGmresItStokes) { |
| 413 |
|
|
// we call the function of the mother class |
| 414 |
|
✗ |
LinearElasticityModel::finalizeForward(); |
| 415 |
|
|
|
| 416 |
|
|
//first we store the new data |
| 417 |
|
✗ |
m_numFixedPointItPerTimeStep.push_back(cIteration); |
| 418 |
|
✗ |
m_totNumOfGMRESItPerTimeStepLE.push_back(cGmresItLE); |
| 419 |
|
✗ |
m_totNumOfGMRESItPerTimeStepStokes.push_back(cGmresItStokes); |
| 420 |
|
|
|
| 421 |
|
✗ |
m_cumulatedForwardTime.push_back(m_timeChecker->diff_cumul()); |
| 422 |
|
✗ |
m_cumulatedStokesTime.push_back(m_t1->diff_cumul()); |
| 423 |
|
✗ |
m_cumulatedLETime.push_back(m_t2->diff_cumul()); |
| 424 |
|
|
|
| 425 |
|
✗ |
writeTimeInfo(); |
| 426 |
|
|
} |
| 427 |
|
|
|
| 428 |
|
|
// we reset the interface quantities of both problem |
| 429 |
|
|
// typically we are going to use this function with |
| 430 |
|
|
// flag = only-rhs |
| 431 |
|
|
// but there is no default |
| 432 |
|
|
void |
| 433 |
|
✗ |
StokesLinearElasticityCoupledModel::clearMatrixRHSOfPbs( FlagMatrixRHS flag) { |
| 434 |
|
✗ |
m_lpbLE->clearMatrixRHS(flag); |
| 435 |
|
✗ |
m_lpbStokes->clearMatrixRHS(flag); |
| 436 |
|
|
} |
| 437 |
|
|
|
| 438 |
|
|
// this function solve the Stokes problem |
| 439 |
|
|
// we do not have a similar function for NS since it is in the mother class |
| 440 |
|
|
void |
| 441 |
|
✗ |
StokesLinearElasticityCoupledModel::solveStokes(int nfp) { |
| 442 |
|
✗ |
if ( FelisceParam::instance().computeSteklov ) { |
| 443 |
|
✗ |
if ( FelisceParam::verbose() ) |
| 444 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Stokes problem using Steklov operator\n"); |
| 445 |
|
✗ |
m_lpbStokes->solveStokesUsingSteklov(m_fstransient, nfp); |
| 446 |
|
|
} else { |
| 447 |
|
✗ |
if ( FelisceParam::verbose() ) |
| 448 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Solving Stokes equations \n"); |
| 449 |
|
|
// In the first fixed point iteration, we assemble both matrix and rhs, consider that for stokes the matrix does not really change since we have no advection. |
| 450 |
|
|
// In the following one we only recompute the the rhs. |
| 451 |
|
|
|
| 452 |
|
|
// Some things to think about it.. |
| 453 |
|
|
// (1) Stokes does not have advection, so the matrix is all in matrix(1) |
| 454 |
|
|
// (2) We cannot have time-depending force term because the reduced steklov method does not allow it. |
| 455 |
|
|
// (3) the only things that can vary with time are the not-constant boundary conditions which are below. |
| 456 |
|
|
|
| 457 |
|
✗ |
if ( nfp == 1 && m_fstransient->iteration == 1) { |
| 458 |
|
✗ |
m_lpbStokes->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::matrix_and_rhs); // important for supg |
| 459 |
|
✗ |
m_lpbStokes->addMatrixRHS(); |
| 460 |
|
✗ |
m_lpbStokes->createAndCopyMatrixRHSWithoutBC(); |
| 461 |
|
✗ |
m_lpbStokes->finalizeEssBCTransient(); // I think it is in preparation of applyBC why then it is not into applyBC? |
| 462 |
|
✗ |
m_lpbStokes->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::matrix_and_rhs,FlagMatrixRHS::matrix_and_rhs); // there is a call here to assemblyNaturalBoundaryConditions. |
| 463 |
|
|
} |
| 464 |
|
✗ |
m_lpbStokes->finalizeEssBCTransient(); // I think it is in preparation of applyBC why then it is not into applyBC? |
| 465 |
|
✗ |
m_lpbStokes->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::only_rhs,FlagMatrixRHS::only_rhs); // there is a call here to assemblyNaturalBoundaryConditions. |
| 466 |
|
|
|
| 467 |
|
✗ |
m_lpbStokes->chronoRS()->start(); |
| 468 |
|
|
|
| 469 |
|
✗ |
m_lpbStokes->solve(MpiInfo::rankProc(), MpiInfo::numProc()); |
| 470 |
|
✗ |
m_lpbStokes->chronoRS()->stop(); |
| 471 |
|
✗ |
m_lpbStokes->exportOutputSteklov(m_fstransient,nfp); |
| 472 |
|
|
} |
| 473 |
|
|
} |
| 474 |
|
|
|
| 475 |
|
|
// two nice banners for the sub-iterations |
| 476 |
|
|
void |
| 477 |
|
✗ |
StokesLinearElasticityCoupledModel::displayBannerIterationInit( felInt cIt ) const { |
| 478 |
|
|
|
| 479 |
|
✗ |
std::stringstream banner; |
| 480 |
|
✗ |
for( felInt k(0); k<25; k++) |
| 481 |
|
✗ |
banner<<"="; |
| 482 |
|
✗ |
banner<<"> iteration "<<cIt<<std::endl; |
| 483 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"%s", banner.str().c_str()); |
| 484 |
|
|
} |
| 485 |
|
|
void |
| 486 |
|
✗ |
StokesLinearElasticityCoupledModel::displayBannerIterationEnd( felInt /*cIt*/, double testQuantity, double testVelocity, double testStresses ) const { |
| 487 |
|
|
|
| 488 |
|
✗ |
std::stringstream banner; |
| 489 |
|
✗ |
for( felInt k(0); k<25; k++) |
| 490 |
|
✗ |
banner<<"="; |
| 491 |
|
✗ |
banner<<"> Test Quantity :"<<testQuantity<<std::endl; |
| 492 |
|
✗ |
for( felInt k(0); k<25; k++) |
| 493 |
|
✗ |
banner<<" "; |
| 494 |
|
✗ |
banner<<" Velocity :"<<testVelocity<<std::endl; |
| 495 |
|
✗ |
for( felInt k(0); k<25; k++) |
| 496 |
|
✗ |
banner<<" "; |
| 497 |
|
✗ |
banner<<" Stresses :"<<testStresses<<std::endl; |
| 498 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"%s", banner.str().c_str()); |
| 499 |
|
|
} |
| 500 |
|
|
} |
| 501 |
|
|
|