GCC Code Coverage Report

Directory: ./
File: Solver/bdf.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 10 20 50.0%
Branches: 0 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 #ifndef _BDF_HPP
16 #define _BDF_HPP
17
18 // System includes
19
20 // External includes
21 #include "Core/NoThirdPartyWarning/Petsc/vec.hpp"
22
23 // Project includes
24 #include "Core/felisce.hpp"
25 #include "PETScInterface/petscVector.hpp"
26
27 namespace felisce
28 {
29 const int BDF_MAX_ORDER = 3;
30
31 /*!
32 \class Bdf
33 \brief Backward differencing formula time discretization
34
35 A differential equation of the form
36
37 \f$ M u' = A u + f \f$
38
39 is discretized in time as
40
41 \f$ M p'(t_{k+1}) = A u_{k+1} + f_{k+1} \f$
42
43 where p denotes the polynomial of order n in t that interpolates
44 (t_i,u_i) for i = k-n+1,...,k+1.
45
46 The approximative time derivative \f$ p'(t_{k+1}) \f$ is a linear
47 combination of state vectors u_i:
48
49 \f$ p'(t_{k+1}) = \frac{1}{\Delta t} (\alpha_0 u_{k+1} - \sum_{i=1}^n \alpha_i u_{k+1-i} )\f$
50
51 Thus we have
52
53 \f$ \frac{\alpha_0}{\Delta t} M u_{k+1} = A u_{k+1} + f + M \bar{p} \f$
54
55 with
56
57 \f$ \bar{p} = \frac{1}{\Delta t} \sum_{i=1}^n \alpha_i u_{k+1-i} \f$
58
59 This class stores the n last state vectors in order to be able to
60 calculate \f$ \bar{p} \f$. It also provides alpha_i
61 and can extrapolate the new state from the n last states with a
62 polynomial of order n-1:
63
64 \f$ u_{k+1} \approx \sum_{i=0}^{n-1} \beta_i u_{k-i} \f$
65 */
66
67 /*!
68 \class Bdf
69 \authors A. Collin
70 \date 12/05/2011
71 \brief ???
72 */
73
74 class Bdf {
75 public:
76 Bdf();
77 ~Bdf();
78 void defineOrder(int n, int nComp=1);
79
80 //!Bdf1.
81 void initialize(const PetscVector& sol_0);
82 //!Bdf2.
83 void initialize(const PetscVector& sol_0,const PetscVector& sol_1);
84 //!Bdf3.
85 void initialize(const PetscVector& sol_0,const PetscVector& sol_1,const PetscVector& sol_2);
86
87 // // //!Bdf1.
88 void initialize(std::vector<PetscVector>& sol_0);
89 // // //!Bdf2.
90 void initialize(std::vector<PetscVector>& sol_0,std::vector<PetscVector>& sol_1);
91 // // //!Bdf3.
92 void initialize(std::vector<PetscVector>& sol_0,std::vector<PetscVector>& sol_1,std::vector<PetscVector>& sol_2);
93
94 void reinitialize(int order, const PetscVector& sol0, const PetscVector& sol1, const PetscVector& sol2);
95
96 void update(PetscVector& sol_n);
97 void update(std::vector<PetscVector>& sol_n);
98
99
100 void computeRHSTime(double dt, PetscVector& RHSTime);
101 void computeRHSTime(double dt, std::vector<PetscVector>& RHSTime);
102 void computeRHSTimeByPart(double dt, std::vector<PetscVector>& RHSTimeByPart);
103 void computeRHSTime(double dt);
104 void extrapolate( PetscVector& extrap);
105 void extrapolate( std::vector<PetscVector>& extrap);
106 void extrapolateByPart(std::vector<PetscVector>& partExtrap);
107
108 //return \alpha_0.
109 inline const double & coeffDeriv0() const {
110 return m_alpha[0];
111 }
112 9592306 inline double & coeffDeriv0() {
113 9592306 return m_alpha[0];
114 }
115
116 // return alpha
117 inline const std::vector<double> & alpha() const {
118 return m_alpha;
119 }
120 inline std::vector<double> & alpha() {
121 return m_alpha;
122 }
123
124 // return beta
125 inline const std::vector<double> & beta() const {
126 return m_beta;
127 }
128 inline std::vector<double> & beta() {
129 return m_beta;
130 }
131
132 // return alpha[i]
133 inline const double & alpha(int i) const {
134 return m_alpha[i];
135 }
136 ✗ inline double & alpha(int i) {
137 ✗ return m_alpha[i];
138 }
139
140 // return beta[i]
141 inline const double & beta(int i) const {
142 return m_beta[i];
143 }
144 ✗ inline double & beta(int i) {
145 ✗ return m_beta[i];
146 }
147
148 //!Access functions.
149
150 inline const std::vector<PetscVector> & vec_sol_n() const {
151 return m_sol_n;
152 }
153 inline std::vector<PetscVector> & vec_sol_n() {
154 return m_sol_n;
155 }
156
157 inline const std::vector<PetscVector> & vec_sol_n_1() const {
158 return m_sol_n_1;
159 }
160 ✗ inline std::vector<PetscVector> & vec_sol_n_1() {
161 ✗ return m_sol_n_1;
162 }
163
164 inline const std::vector<PetscVector> & vec_sol_n_2() const {
165 return m_sol_n_2;
166 }
167 ✗ inline std::vector<PetscVector> & vec_sol_n_2() {
168 ✗ return m_sol_n_2;
169 }
170
171 inline const std::vector<PetscVector> & vec_RHS() const {
172 return m_rhs;
173 }
174 inline std::vector<PetscVector> & vec_RHS() {
175 return m_rhs;
176 }
177
178 inline const PetscVector& sol_n() const {
179 return m_sol_n[0];
180 }
181 31 inline PetscVector& sol_n() {
182 31 return m_sol_n[0];
183 }
184
185 inline const PetscVector& sol_n_1() const {
186 return m_sol_n_1[0];
187 }
188 20 inline PetscVector& sol_n_1() {
189 20 return m_sol_n_1[0];
190 }
191
192 inline const PetscVector& sol_n_2() const {
193 return m_sol_n_2[0];
194 }
195 ✗ inline PetscVector& sol_n_2() {
196 ✗ return m_sol_n_2[0];
197 }
198
199 inline const PetscVector& vector() const {
200 return m_rhs[0];
201 }
202 2156 inline PetscVector& vector() {
203 2156 return m_rhs[0];
204 }
205
206 inline const int & order() const {
207 return m_order;
208 }
209 140 inline int & order() {
210 140 return m_order;
211 }
212
213 inline const int & numComp() const {
214 return m_numberOfComp;
215 }
216 inline int & numComp() {
217 return m_numberOfComp;
218 }
219
220 // Access fonction to use it in user file
221 inline const std::vector<PetscVector> & vec_solExt() const {
222 return m_solExtrapol;
223 }
224 inline std::vector<PetscVector> & vec_solExt() {
225 return m_solExtrapol;
226 }
227 inline const PetscVector& solExt() const {
228 return m_solExtrapol[0];
229 }
230 inline PetscVector& solExt() {
231 return m_solExtrapol[0];
232 }
233
234 //! Returns the time derivative of the solution with u^{n+1}=vec
235 void time_der( double dt, const PetscVector& m_vecSol, PetscVector& deriv ) const;
236 void time_der( double dt, const std::vector<PetscVector>& m_vecSol, std::vector<PetscVector>& deriv ) const;
237
238
239 private:
240 //! Order of the BDF derivative/extrapolation: the time-derivative
241 //! coefficients std::vector has size n+1, the extrapolation std::vector has size n
242 int m_order;
243 //! Number of components
244 int m_numberOfComp;
245 //! Size du vecteur solution
246 felInt m_size;
247 //! Coefficients \f$ \alpha_i \f$ of the time bdf discretization
248 std::vector<double> m_alpha;
249 //! Coefficients \f$ \beta_i \f$ of the extrapolation
250 std::vector<double> m_beta;
251
252 std::vector<PetscVector> m_sol_n;
253 std::vector<PetscVector> m_sol_n_1;
254 std::vector<PetscVector> m_sol_n_2;
255 std::vector<PetscVector> m_rhs;
256
257 // extrapolated solution when using bdf order > 1
258 std::vector<PetscVector> m_solExtrapol;
259
260 bool m_build_order_2;
261 bool m_build_order_3;
262 bool m_build_RHS;
263 };
264 }
265
266 #endif
267