GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemPoroElasticity.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 127 300 42.3%
Branches: 135 614 22.0%
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/linearProblemPoroElasticity.hpp"
21 #include "FiniteElement/elementVector.hpp"
22 #include "FiniteElement/elementMatrix.hpp"
23 #include "InputOutput/io.hpp"
24
25 namespace felisce {
26 2 LinearProblemPoroElasticity::LinearProblemPoroElasticity():
27
7/14
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
 2 LinearProblem("Linear Poro Elasticity",/*number of matrices*/2) {
28
29 // Default value for the initial (constant) pressure
30 2 m_pInit=0;
31
32 //========Initialization of the PetscVectors=============//
33 2 std::vector<std::string> listPar,listSeq;
34
35
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listPar.emplace_back("solutionPreviousTimeStep");
36
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listSeq.emplace_back("solutionPreviousTimeStep");
37
38
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listPar.emplace_back("filtrationVelocity");
39
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listSeq.emplace_back("filtrationVelocity");
40
41 // We do not need the sequential counterpart.
42
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listPar.emplace_back("massTimesPressureGradient");
43
44
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listPar.emplace_back("pressureGradient");
45
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listSeq.emplace_back("pressureGradient");
46
47 // We do not need the sequential counterpart;
48
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listPar.emplace_back("anisotropicDirection");
49
50
2/2
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
 12 for( std::size_t k(0); k<listPar.size();++k) {
51
2/4
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
 10 m_vecs.Init(listPar[k]);
52 }
53
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
 8 for( std::size_t k(0); k<listSeq.size();++k) {
54
2/4
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
 6 m_seqVecs.Init(listSeq[k]);
55 }
56
57
58 // Default value
59 2 m_matricesForFiltrationVelocityBuilt=false;
60
61 //==============Scaling the value of the boundary conditions========//
62 // Not tested too much.
63
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
 2 if ( FelisceParam::instance().scalePressureEq ) {
64 ✗ int c(0);
65 ✗ for ( std::size_t k(0); k<FelisceParam::instance().variable.size(); ++k) {
66 ✗ std::cout<<"Considering BC number: "<<k
67 ✗ <<" var: "<<FelisceParam::instance().variable[k]
68 ✗ <<" comp: "<<FelisceParam::instance().component[k]
69 ✗ <<" type: "<<FelisceParam::instance().type[k]
70 ✗ <<" value: "<<FelisceParam::instance().value[k]
71 ✗ <<std::endl;
72 ✗ int n(2);
73 ✗ if ( FelisceParam::instance().component[k] == "Comp123" )
74 ✗ n=3;
75 ✗ if ( FelisceParam::instance().component[k] == "Comp1" ||
76 ✗ FelisceParam::instance().component[k] == "CompNA" )
77 ✗ n=1;
78
79 ✗ if ( FelisceParam::instance().variable[k] == "pressure" ) {
80 ✗ if ( FelisceParam::instance().type[k] == "Neumann" ) {
81 ✗ std::cout<<"Rescaling "<<k<<"th bc. Values from "<<c<<" to "<<c+n-1<<std::endl;
82 ✗ for ( int j(0); j<n; ++j ) {
83 ✗ std::cout<<FelisceParam::instance().value[c+j]<<" --> ";
84 ✗ FelisceParam::instance().value[c+j] *= FelisceParam::instance().timeStep*FelisceParam::instance().biotModulus;
85 ✗ std::cout<<FelisceParam::instance().value[c+j]<<std::endl;
86 }
87 } else {
88 ✗ std::cout<<"Bc n."<<k<<" was not rescaled because it was not of type neumann"<<std::endl;
89 }
90 }
91 ✗ c=c+n;
92 }
93 }
94 2 }
95 2 void LinearProblemPoroElasticity::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) {
96
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 LinearProblem::initialize(mesh,comm,doUseSNES);
97
98 //=================Setting the physical parameters=================//
99 2 m_fstransient=fstransient;
100
101
5/10
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
 2 if ( FelisceParam::instance().young < 0 || FelisceParam::instance().poisson < 0 ) {
102 ✗ m_mu = FelisceParam::instance().mu_lame;
103 ✗ m_lambda = FelisceParam::instance().lambda_lame;
104 } else {
105
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 m_mu = FelisceParam::instance().young/2/(1+FelisceParam::instance().poisson);
106
4/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
 2 m_lambda = FelisceParam::instance().young*FelisceParam::instance().poisson/(1.+FelisceParam::instance().poisson)/(1.-2.*FelisceParam::instance().poisson);
107 }
108
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_deltaT = FelisceParam::instance().timeStep;
109
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_b = FelisceParam::instance().biotCoefficient;
110
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_M = FelisceParam::instance().biotModulus;
111
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_k0 = FelisceParam::instance().isoPermeab;
112
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_k1 = FelisceParam::instance().anisoPermeab;
113
114
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 this->displayPoroElasticProperties();
115
116
117 //================Initializing variables and unknowns=============//
118 2 std::vector<PhysicalVariable> listVariable;
119 2 std::vector<std::size_t> listNumComp;
120
121 // Displacement
122
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listVariable.push_back(displacement);
123
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 listNumComp.push_back(this->dimension());
124 // Pore-pressure
125
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listVariable.push_back(pressure);
126
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listNumComp.push_back(1);
127
128
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 userAddOtherVariables(listVariable,listNumComp);
129
130 // Declare them as unknown of the linear system
131
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_listUnknown.push_back(displacement);
132
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_listUnknown.push_back(pressure);
133
134 // Declare them as physical variable
135
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 this->definePhysicalVariable(listVariable,listNumComp);
136
137 // Save the variable id of the displacement and of the pressure
138
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 m_iDisplacement = this->listVariable().getVariableIdList(displacement);
139
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 m_iPressure = this->listVariable().getVariableIdList(pressure);
140 2 }
141
142 2 void LinearProblemPoroElasticity::displayPoroElasticProperties() const {
143 // drained parameters
144 2 double shearD = m_mu;
145 2 double lambdaD = m_lambda;
146 2 double bulkD = lambdaD+2./3.*shearD;
147 2 double poissonD = (3*bulkD-2*shearD)/2/(3*bulkD+shearD);
148 2 double youngD = 9*bulkD*shearD/(3*bulkD+shearD);
149
150 // undrained parameters
151 2 double shearU = shearD;
152 2 double lambdaU = lambdaD + m_b*m_b*m_M;
153 2 double bulkU = bulkD + m_b*m_b*m_M;
154 2 double poissonU = lambdaU / (2*(lambdaU+shearD));
155 2 double youngU = (1+poissonU)/(1+poissonD)*youngD;
156
157
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 std::stringstream msg;
158
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 msg<<"========================================" << std::endl;
159
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 msg<<"--- DRAINED PROPERTIES (p=0):" << std::endl;
160
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Shear Modulus (G or mu): "<< shearD << std::endl;
161
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Lame's coefficient (lambda): "<< lambdaD << std::endl;
162
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Bulk modulus (K): "<< bulkD << std::endl;
163
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Young Modulus: " << youngD << std::endl;
164
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Poisson ratio: "<< poissonD << std::endl;
165
166
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 msg<<"--- UNDRAINED PROPERTIES (no extra mass):" << std::endl;
167
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Shear Modulus (G or mu): "<< shearU << std::endl;
168
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Lame's coefficient (lambda): "<< lambdaU << std::endl;
169
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Bulk modulus (K): "<< bulkU << std::endl;
170
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Young Modulus: " << youngU << std::endl;
171
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Poisson ratio: "<< poissonU << std::endl;
172
173
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 msg<<"--- OTHER PARAMETERS-----"<<std::endl;
174
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Time step: "<< m_deltaT << std::endl;
175
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Biot coefficient(b): "<< m_b << std::endl;
176
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Biot modulus (M): "<< m_M << std::endl;
177
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Isotropic permeability (k0)"<< m_k0<<std::endl;
178
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 msg<<" Ansotropic permeability (k1)"<< m_k1<<std::endl;
179
180
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 2 PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
181 2 }
182 22 void LinearProblemPoroElasticity::initPerElementType(ElementType /*eltType*/, FlagMatrixRHS /*flagMatrixRHS*/) {
183 // Pointers to the current finite element.
184 22 m_feDisp = m_listCurrentFiniteElement[m_iDisplacement];
185 22 m_fePre = m_listCurrentFiniteElement[m_iPressure];
186
187 // Element fields to store the old values of displacement and pore-pressure
188 22 m_elemFieldDisplacement.initialize(DOF_FIELD,*m_feDisp,this->dimension());
189 22 m_elemFieldPressure.initialize(DOF_FIELD,*m_fePre,1);
190 22 }
191
192 18 void LinearProblemPoroElasticity::advanceInTime() {
193
194 // User dependent computations.
195 18 userPostProcessing();
196
197
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
 18 if ( FelisceParam::instance().exportFiltrationVelocity) {
198 // Computation of filtration velocity.
199 ✗ if ( ! m_matricesForFiltrationVelocityBuilt ) {
200 ✗ this->assembleMatricesForFiltrationVelocity();
201 ✗ m_matricesForFiltrationVelocityBuilt = true;
202 }
203 ✗ this->computeFiltrationVelocity();
204 }
205 // Move forward
206
3/6
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 18 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 18 times.
✗ Branch 10 not taken.
 18 m_vecs.Get("solutionPreviousTimeStep").copyFrom(this->solution());
207
2/4
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 18 times.
✗ Branch 6 not taken.
 18 this->gather("solutionPreviousTimeStep");
208 18 }
209
210 2662 void LinearProblemPoroElasticity::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
211
212 // UPDATE PHASE:
213 // -- of the finite element
214 2662 m_feDisp->updateFirstDeriv(0, elemPoint);
215 2662 m_fePre->updateFirstDeriv(0, elemPoint);
216 // -- of the anisotropic permeability direction
217 2662 updateAnisotropicDirection(elemPoint);
218 // -- of the elementfields
219
2/2
✓ Branch 1 taken 242 times.
✓ Branch 2 taken 2420 times.
 2662 if (m_fstransient->iteration == 1) { // Setting intial condition on the pressure
220 242 m_elemFieldPressure.setValue(m_pInit,0);
221 } else {
222
3/6
✓ Branch 3 taken 2420 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2420 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2420 times.
✗ Branch 10 not taken.
 2420 m_elemFieldPressure.setValue(m_seqVecs.Get("solutionPreviousTimeStep") ,*m_fePre,iel,m_iPressure,m_ao,dof());
223 }
224
3/6
✓ Branch 3 taken 2662 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2662 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2662 times.
✗ Branch 10 not taken.
 2662 m_elemFieldDisplacement.setValue(m_seqVecs.Get("solutionPreviousTimeStep"),*m_feDisp,iel,m_iDisplacement,m_ao,dof());
225
226
227 // Wether the pressure equation should be multiplied by (dt*M) or by one.
228 2662 double scalingCoeff(1.);
229
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2662 times.
 2662 if ( FelisceParam::instance().scalePressureEq )
230 ✗ scalingCoeff = m_deltaT*m_M;
231
232 2662 felInt pressureBlock( this->dimension() );
233 2662 felInt displacementBlock(0);
234
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2662 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2662 if ( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix)) {
235
2/2
✓ Branch 1 taken 242 times.
✓ Branch 2 taken 2420 times.
 2662 if ( m_fstransient->iteration == 0 ) {
236 242 std::size_t idMat = 1;// Everything should be placed in the static matrix (idMat=1)
237 //! ==================================DISPLACEMENT EQUATION=======================================================
238 //! Displacement equation: quasi-static, no mass here!
239 //! Divergence of the elastic component of the stress tensor just as in linear elasticity
240 242 m_elementMat[idMat]->eps_phi_i_eps_phi_j( 2*m_mu, *m_feDisp, /*iblock*/displacementBlock, /*jblock*/displacementBlock, /*numOfBlockRepetitions*/this->dimension());
241 242 m_elementMat[idMat]->div_phi_j_div_phi_i( m_lambda, *m_feDisp, /*iblock*/displacementBlock, /*jblock*/displacementBlock, /*numOfBlockRepetitions*/this->dimension());
242 //! now the pressure part
243 242 m_elementMat[idMat]->psi_j_div_phi_i(-m_b,*m_feDisp,*m_fePre , /*iblock*/displacementBlock, /*jblock*/pressureBlock);
244
245 //! ==================================PRESSURE EQUATION===========================================================
246 //! ------- First term coming from the time derivative of the divergence of the displacement
247 242 m_elementMat[idMat]->psi_i_div_phi_j(scalingCoeff*m_b/m_deltaT,*m_feDisp,*m_fePre, /*iblock*/pressureBlock,/*jblock*/displacementBlock);
248 //! ------- First term coming from the time derivative of the pressure
249 242 m_elementMat[idMat]->phi_i_phi_j(scalingCoeff/(m_M*m_deltaT),*m_fePre,/*iblock*/pressureBlock,/*jblock*/pressureBlock,/*numOfBlockRepetitions*/1);
250
251 // Pressure equation: laplacian terms.
252 // Anisotropic part of the diffusive term
253 242 m_elementMat[idMat]->tau_grad_phi_i_tau_grad_phi_j(scalingCoeff*m_k1,m_anisotropicDir,*m_fePre,/*iblock*/pressureBlock,/*jblock*/pressureBlock,/*nbOfBlockRep*/1,/*domain dim*/this->dimension());
254 // Isotropic part
255 242 m_elementMat[idMat]->grad_phi_i_grad_phi_j(scalingCoeff*m_k0,*m_fePre,/*iblock*/pressureBlock, /*jblock*/pressureBlock,/*numOfBlockRepetitions*/1);
256 }
257 }
258
259 // Right hand side:
260
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2662 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2662 if ( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs)) {
261 // Pressure equation:
262 // ------- Second term coming from the time derivative of the divergence of the displacement
263 2662 m_elementVector[0]->div_f_phi_i(scalingCoeff*m_b/m_deltaT, *m_feDisp,*m_fePre,m_elemFieldDisplacement,/*iblock*/pressureBlock);
264 // ------- Second term coming from the time derivative of the pressure
265 2662 m_elementVector[0]->source(scalingCoeff/(m_M*m_deltaT), *m_fePre,m_elemFieldPressure,/*iblock*/pressureBlock,/*numOfBlockRepetitions*/1);
266 }
267 // OTHER TERMS:
268 // Mass source and external load should be coded in userComputeElementArray
269 // Boundary terms such as Neumann force terms should be standard and therefore implemented in linearProblem.cpp
270 // Pay attention:
271 // if the pressure equation is multiplied by m_M you have to be consistent when coding the user files.
272 2662 }
273 2662 void LinearProblemPoroElasticity::updateAnisotropicDirection(const std::vector<Point*>& /*elemPoint*/) {
274 // By default there is no anisotropic permeability direction.
275 // --> Customize this function in the user
276
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2660 times.
 2662 if ( m_anisotropicDir.size() == 0) {
277 2 m_anisotropicDir.resize(this->dimension()*m_fePre->numDof());
278 } else {
279
2/2
✓ Branch 1 taken 15960 times.
✓ Branch 2 taken 2660 times.
 18620 for (std::size_t i(0); i < m_anisotropicDir.size(); ++i) {
280 15960 m_anisotropicDir[i]=0.0;
281 }
282 }
283 2662 }
284 // Exporting a zero-initial filtration velocity.
285 2 void LinearProblemPoroElasticity::exportInitialFiltrationVelocity(std::vector<IO::Pointer>& io) {
286
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if(FelisceParam::instance().exportFiltrationVelocity) {
287 ✗ PetscVector zero;
288 ✗ zero.duplicateFrom(this->solution());
289 ✗ zero.zeroEntries();
290 ✗ std::stringstream msg;
291 ✗ msg<<"Writing filtration velocity"<<std::endl;
292 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
293 ✗ this->writeSolutionFromVec(zero, MpiInfo::rankProc(), io, m_fstransient->time, m_fstransient->iteration, std::string("filtrationVelocity"));
294 // this->writeSolutionFromVec( zero, MpiInfo::rankProc(), io, m_fstransient->time, m_fstransient->iteration, std::string("anisotropicDirection"));
295 }
296 2 }
297
298 18 void LinearProblemPoroElasticity::exportFiltrationVelocity(std::vector<IO::Pointer>& io) {
299
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
 18 if( FelisceParam::instance().exportFiltrationVelocity ) {
300 ✗ std::stringstream msg;
301 ✗ msg<<"Writing filtration velocity"<<std::endl;
302 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
303 ✗ this->writeSolutionFromVec(m_seqVecs.Get("filtrationVelocity"),
304 ✗ MpiInfo::rankProc(),
305 io,
306 ✗ m_fstransient->time,
307 ✗ m_fstransient->iteration,
308 ✗ std::string("filtrationVelocity"));
309 // this->writeSolutionFromVec(m_vecs.Get("anisotropicDirection"),
310 // MpiInfo::rankProc(),
311 // io,
312 // m_fstransient->time,
313 // m_fstransient->iteration,
314 // std::string("anisotropicDirection"));
315 }
316 18 }
317
318 ✗ void LinearProblemPoroElasticity::assembleMatricesForFiltrationVelocity() {
319 // Vector/Matrices intialization
320 ✗ m_mass.duplicateFrom(matrix(0),MAT_DO_NOT_COPY_VALUES);
321 ✗ m_mass.zeroEntries();
322 ✗ m_derivative.duplicateFrom(matrix(0),MAT_DO_NOT_COPY_VALUES);
323 ✗ m_derivative.zeroEntries();
324
325 ✗ const auto& r_instance = FelisceParam::instance(this->instanceIndex());
326 ✗ const std::string solver = r_instance.solver[m_identifier_solver];
327 ✗ const std::string preconditioner = r_instance.preconditioner[m_identifier_solver];
328 ✗ const double relativeTolerance = r_instance.relativeTolerance[m_identifier_solver];
329 ✗ const double absoluteTolerance = r_instance.absoluteTolerance[m_identifier_solver];
330 ✗ const int maxIteration = r_instance.maxIteration[m_identifier_solver];
331 ✗ const int gmresRestart = r_instance.gmresRestart[m_identifier_solver];
332 ✗ const bool initPrevSolution = r_instance.initSolverWithPreviousSolution[m_identifier_solver];
333 ✗ m_kspForFiltrationVelocity->init();
334 ✗ m_kspForFiltrationVelocity->setKSPandPCType(solver, preconditioner);
335 ✗ m_kspForFiltrationVelocity->setTolerances(relativeTolerance, absoluteTolerance, maxIteration, gmresRestart);
336 ✗ m_kspForFiltrationVelocity->setKSPOptions(solver, initPrevSolution);
337
338 // Loop over the domain
339 // use to define a "global" numbering of element in the mesh.
340 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
341 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
342 ✗ ElementType eltType = (ElementType)ityp;
343 ✗ numElement[eltType] = 0;
344 }
345
346 ✗ std::vector<Point*> elemPoint;
347 ✗ std::vector<felInt> elemIdPoint;
348 felInt ielSupportDof;
349
350 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
351 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
352 ✗ ElementType eltType = bagElementTypeDomain[i];
353
354 // resize array.
355 ✗ int numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
356 ✗ elemPoint.resize(numPointPerElt, nullptr);
357 ✗ elemIdPoint.resize(numPointPerElt, 0);
358
359 ✗ defineFiniteElement(eltType);
360 ✗ initElementArray();
361 ✗ allocateArrayForAssembleMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
362
363 // virtual function use in derived problem to allocate elemenField necessary.
364 ✗ initPerElementType(eltType, FlagMatrixRHS::only_matrix);
365 // second loop on region of the mesh.
366 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin();
367 ✗ itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
368 // int currentLabel = itRef->first;
369 ✗ felInt numEltPerLabel = itRef->second.second;
370 // Third loop on element in the region with the type: eltType. ("real" loop on elements).
371 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
372 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, &ielSupportDof);
373
374 // FE update
375 ✗ m_listCurrentFiniteElement[1/*pressure*/]->updateFirstDeriv(0, elemPoint);
376
377 // MASS
378 ✗ m_elementMat[0]->zero();
379 ✗ m_elementMat[0]->phi_i_phi_j(1., *m_listCurrentFiniteElement[1],0,0,this->dimension()+1);
380 ✗ setValueCustomMatrix(ielSupportDof,ielSupportDof,m_mass);
381
382 // Derivative
383 ✗ m_elementMat[0]->zero();
384 ✗ for (int h(0); h<this->dimension();++h) {
385 ✗ m_elementMat[0]->dpsi_j_dh_psi_i(1.0,*m_listCurrentFiniteElement[1],*m_listCurrentFiniteElement[1], h, h, this->dimension());
386 }
387 ✗ setValueCustomMatrix(ielSupportDof,ielSupportDof,m_derivative);
388
389 ✗ numElement[eltType]++;
390 }
391 }
392 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
393 }
394 ✗ m_mass.assembly(MAT_FINAL_ASSEMBLY);
395 ✗ m_derivative.assembly(MAT_FINAL_ASSEMBLY);
396
397 ✗ const std::string preconditionerOption = r_instance.setPreconditionerOption[m_identifier_solver];
398 ✗ m_kspForFiltrationVelocity->setKSPOperator(m_mass, preconditionerOption);
399 }
400
401 ✗ void LinearProblemPoroElasticity::computeFiltrationVelocity() {
402 ✗ std::stringstream msg;
403 ✗ msg<<"Solution of the linear system to compute the filtration velocity: "<<std::endl;
404 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
405 ✗ mult(m_derivative,this->solution(),m_vecs.Get("massTimesPressureGradient")); // aux_i = < \nabla p, \phi_i >
406
407 ✗ const auto& r_instance = FelisceParam::instance(this->instanceIndex());
408 ✗ m_kspForFiltrationVelocity->solve(m_vecs.Get("massTimesPressureGradient"), m_vecs.Get("pressureGradient"), r_instance.linearSolverVerbose);
409 ✗ this->gather("pressureGradient");
410
411 // correctly Perform the product v=Kv;
412 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
413 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
414 ✗ ElementType eltType = (ElementType)ityp;
415 ✗ numElement[eltType] = 0;
416 }
417
418 ✗ std::vector<Point*> elemPoint;
419 ✗ std::vector<felInt> elemIdPoint;
420 felInt ielSupportDof;
421
422 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
423 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
424 ✗ ElementType eltType = bagElementTypeDomain[i];
425
426 // resize array.
427 ✗ int numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
428 ✗ elemPoint.resize(numPointPerElt, nullptr);
429 ✗ elemIdPoint.resize(numPointPerElt, 0);
430
431 ✗ defineFiniteElement(eltType);
432 ✗ initElementArray();
433 ✗ allocateArrayForAssembleMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
434
435 // virtual function use in derived problem to allocate elemenField necessary.
436 ✗ initPerElementType(eltType, FlagMatrixRHS::only_matrix);
437 // second loop on region of the mesh.
438 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin();
439 ✗ itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
440 // int currentLabel = itRef->first;
441 ✗ felInt numEltPerLabel = itRef->second.second;
442 // Third loop on element in the region with the type: eltType. ("real" loop on elements).
443 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
444 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, &ielSupportDof);
445 ✗ this->updateAnisotropicDirection(elemPoint);
446
447 ✗ m_elementVector[0]->zero();
448 ✗ felInt counter(0);//we assume the displacement is the first unknown
449 ✗ for (int icomp(0); icomp<this->dimension();++icomp) {
450 ✗ for ( int idof(0); idof<m_fePre->numDof(); ++idof) {
451 ✗ m_elementVector[0]->vec()(counter) = m_anisotropicDir[this->dimension()*idof+icomp];
452 ✗ counter++;
453 }
454 }
455 ✗ setValueCustomVector(ielSupportDof,INSERT_VALUES,m_vecs.Get("anisotropicDirection"));
456
457 ✗ m_elementVector[0]->zero();
458 ✗ counter=0;//we assume the displacement is the first unknown
459 ✗ for (int icomp(0); icomp<this->dimension();++icomp) {
460 ✗ for ( int idof(0); idof<m_fePre->numDof(); ++idof) {
461 //identify physical point
462 double x,y,z;
463 ✗ if ( (std::size_t)idof >= elemPoint.size() ) { // ONLY FOR P2, assuming prisms
464 ✗ x=0.5*(elemPoint[idof-3]->x()+elemPoint[idof-6]->x());
465 ✗ y=0.5*(elemPoint[idof-3]->y()+elemPoint[idof-6]->y());
466 ✗ z=0.5*(elemPoint[idof-3]->z()+elemPoint[idof-6]->z());
467 } else {
468 ✗ x=elemPoint[idof]->x();
469 ✗ y=elemPoint[idof]->y();
470 ✗ z=elemPoint[idof]->z();
471 }
472 ✗ for (int jcomp(0); jcomp<this->dimension();++jcomp) {
473 //compute -phi(idof)*tensor(idof)_{icomp,jcomp}
474 ✗ double val=m_k1*m_anisotropicDir[this->dimension()*idof+icomp]*m_anisotropicDir[this->dimension()*idof+jcomp];
475 ✗ if (icomp==jcomp)
476 ✗ val+=m_k0;
477 ✗ val = -val*this->userEvaluatePorosity(x,y,z);
478
479 //retrieving the good value from pressureGradient
480 double valOfGradP;
481 felInt gdof;
482 ✗ dof().loc2glob(ielSupportDof,idof,0/*idVariable*/,jcomp,gdof);
483 ✗ AOApplicationToPetsc(ao(),1,&gdof);
484 ✗ m_seqVecs.Get("pressureGradient").getValues(1,&gdof,&valOfGradP);
485
486 //summing
487 ✗ m_elementVector[0]->vec()(counter) += val*valOfGradP;
488 }
489 ✗ counter++;
490 }
491 }
492 ✗ setValueCustomVector(ielSupportDof,INSERT_VALUES,m_vecs.Get("filtrationVelocity"));
493 ✗ numElement[eltType]++;
494 }
495 }
496 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
497 }
498 ✗ m_vecs.Get("filtrationVelocity").assembly();
499 ✗ this->gather("filtrationVelocity");
500 }
501
502 #ifdef FELISCE_WITH_CVGRAPH
503 /*! \brief this function assemble the mass matrix at the boundary.
504 * And the KSP
505 */
506 ✗ void LinearProblemPoroElasticity::assembleMassBoundaryAndInitKSP( std::size_t iBD ) {
507 ✗ LinearProblem::assembleMassBoundaryAndInitKSP( &LinearProblemPoroElasticity::massMatrixComputer,
508 ✗ FelisceParam::instance().interfaceLabels,
509 &LinearProblemPoroElasticity::initPerETMass,
510 &LinearProblemPoroElasticity::updateFE,
511 iBD);
512 }
513
514
515
516 /*! \brief Function to assemble the mass matrix
517 * Function to be called in the assemblyLoopBoundaryGeneral
518 */
519 ✗ void LinearProblemPoroElasticity::massMatrixComputer(felInt ielSupportDof) {
520 ✗ this->m_elementMatBD[0]->zero();
521 ✗ this->m_elementMatBD[0]->phi_i_phi_j(/*coef*/1.,*m_curvFeDisp,/*iblock*/0,/*iblock*/0,this->dimension());
522 ✗ this->setValueMatrixBD(ielSupportDof);
523 }
524 void
525 ✗ LinearProblemPoroElasticity::initPerETMass() {
526 ✗ felInt numDofTotal = 0;
527 //pay attention this numDofTotal is bigger than expected since it counts for all the VARIABLES
528 ✗ for (std::size_t iFe = 0; iFe < this->m_listCurvilinearFiniteElement.size(); iFe++) {//this loop is on variables while it should be on unknown
529 ✗ numDofTotal += this->m_listCurvilinearFiniteElement[iFe]->numDof()*this->m_listVariable[iFe].numComponent();
530 }
531 ✗ m_globPosColumn.resize(numDofTotal); m_globPosRow.resize(numDofTotal); m_matrixValues.resize(numDofTotal*numDofTotal);
532
533 ✗ m_curvFeDisp = this->m_listCurvilinearFiniteElement[ m_iDisplacement ];
534 }
535 void
536 ✗ LinearProblemPoroElasticity::updateFE(const std::vector<Point*>& elemPoint, const std::vector<int>&) {
537 ✗ m_curvFeDisp->updateMeasNormal(0, elemPoint);
538 }
539 #endif
540 }
541