LifeV
StabilizationSUPGALE.cpp
Go to the documentation of this file.
1 #include <lifev/navier_stokes_blocks/solver/StabilizationSUPGALE.hpp>
2 
3 // MACRO TO DEFINE TAU_M
4 #define TAU_M value(1)/( eval(squareroot,TAU_M_DEN) )
5 #define TAU_M_DEN TAU_M_DEN_DT + TAU_M_DEN_VEL + TAU_M_DEN_VISC
6 #define TAU_M_DEN_DT value(M_density*M_density)*value(M_bdfOrder*M_bdfOrder)/value(M_timestep * M_timestep)
7 #define TAU_M_DEN_VEL value(M_density*M_density)*dot(value(M_fespaceUETA, beta_km1), G * value(M_fespaceUETA, beta_km1))
8 #define TAU_M_DEN_VISC ( value(M_C_I)*value(M_viscosity*M_viscosity) *dot (G, G) )
9 
10 // MACRO TO DEFINE TAU_C
11 #define TAU_C value(1.0)/( dot(g, TAU_M*g ) )
12 
13 namespace LifeV
14 {
15 
16 //=============================================================================
17 // Constructor
18 //=============================================================================
19 
20 StabilizationSUPGALE::StabilizationSUPGALE():
21  M_label("SUPG-ALE")
22 {
23 }
24 
25 //=============================================================================
26 // Methods
27 //=============================================================================
28 
29 void StabilizationSUPGALE::setConstant(const int & value)
30 {
31  if ( value == 1 )
32  M_C_I = 30;
33  else if ( value == 2 )
34  M_C_I = 60;
35  else
36  ASSERT(0!=0, "Please implement a suitable value for M_C_I for your velocity FE order");
37 }
38 
39 void StabilizationSUPGALE::buildGraphs()
40 {
41  // We assume that the fluid velocity and the mesh velocity are discretized by FE of same order
42  vector_Type u_km1( M_uFESpace->map(), Repeated);
43  vector_Type beta_km1( M_uFESpace->map(), Repeated);
44  vector_Type p_km1( M_pFESpace->map(), Repeated);
45  vector_Type u_bdf( M_uFESpace->map(), Repeated);
46 
47  u_km1.zero();
48  beta_km1.zero();
49  p_km1.zero();
50  u_bdf.zero();
51 
52  u_km1 += 1;
53  beta_km1 += 1;
54  p_km1 += 1;
55  u_bdf += 1;
56 
57  std::shared_ptr<SquareRoot_SUPGALE> squareroot(new SquareRoot_SUPGALE());
58 
59  MatrixSmall<3, 3> Eye;
60  Eye *= 0.0;
61  Eye[0][0] = 1;
62  Eye[1][1] = 1;
63  Eye[2][2] = 1;
64 
65  {
66  using namespace ExpressionAssembly;
67 
68  // Graph for block 00
69  M_graph_block00.reset (new Epetra_FECrsGraph (Copy, * (M_uFESpace->map().map (Unique) ), 0) );
70  buildGraph ( elements ( M_uFESpace->mesh() ),
71  quadRuleTetra4pt,
72  M_fespaceUETA,
73  M_fespaceUETA,
74  /*(1)*/ TAU_M * value(M_density*M_density) * value(M_alpha/M_timestep) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), phi_j )
75  /*(2)*/ +TAU_M * value(M_density*M_density) * value(M_alpha/M_timestep) * dot( phi_j*grad(phi_i), value(M_fespaceUETA, u_km1) )
76  /*(3)*/ -TAU_M * value(M_density*M_density) * dot( phi_j*grad(phi_i), value(M_fespaceUETA, u_bdf) )
77  /*(4)*/ +TAU_M * value(M_density*M_density) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), phi_j*grad(M_fespaceUETA, u_km1) )
78  /*(5)*/ +TAU_M * value(M_density*M_density) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), value(M_fespaceUETA, beta_km1)*grad(phi_j) )
79  /*(6)*/ +TAU_M * value(M_density*M_density) * dot( phi_j*grad(phi_i), value(M_fespaceUETA, beta_km1)*grad(M_fespaceUETA, u_km1) )
80  /*(8)*/ +TAU_M * value(M_density) * dot( phi_j*grad(phi_i), grad(M_fespacePETA, p_km1) )
81  /*(9)*/ -TAU_M * value(M_density*M_viscosity) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), laplacian(phi_j) )
82  /*(10)*/ -TAU_M * value(M_density*M_viscosity) * dot( phi_j*grad(phi_i), laplacian(M_fespaceUETA, u_km1) )
83  /*(16)*/ +TAU_C * div(phi_i)*div(phi_j)
84  ) >> M_graph_block00;
85  M_graph_block00->GlobalAssemble();
86  M_graph_block00->OptimizeStorage();
87 
88  // Graph for block 10
89  M_graph_block10.reset (new Epetra_FECrsGraph (Copy, * (M_pFESpace->map().map (Unique) ), 0) );
90  buildGraph ( elements (M_pFESpace->mesh() ),
91  quadRuleTetra4pt,
92  M_fespacePETA,
93  M_fespaceUETA,
94  /*(11)*/ TAU_M * value(M_alpha*M_density/M_timestep) * dot( grad(phi_i), phi_j )
95  /*(12)*/ +TAU_M * value(M_density) * dot( grad(M_fespaceUETA, u_km1)*grad(phi_i), phi_j )
96  /*(13)*/ +TAU_M * value(M_density) * dot( grad(phi_i), value(M_fespaceUETA, beta_km1)*grad(phi_j) )
97  /*(15)*/ -TAU_M * value(M_viscosity) * dot( grad(phi_i), laplacian(phi_j) )
98  ) >> M_graph_block10;
99  M_graph_block10->GlobalAssemble ( * (M_uFESpace->map().map (Unique) ), * (M_pFESpace->map().map (Unique) ) );
100  M_graph_block10->OptimizeStorage();
101 
102  // Graph for block 01
103  M_graph_block01.reset (new Epetra_FECrsGraph (Copy, * (M_uFESpace->map().map (Unique) ), 0) );
104  buildGraph ( elements (M_pFESpace->mesh() ),
105  quadRuleTetra4pt,
106  M_fespaceUETA,
107  M_fespacePETA,
108  /*(7)*/ TAU_M * value(M_density) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), grad(phi_j) )
109  ) >> M_graph_block01;
110  M_graph_block01->GlobalAssemble ( * (M_pFESpace->map().map (Unique) ), * (M_uFESpace->map().map (Unique) ) );
111  M_graph_block01->OptimizeStorage();
112 
113  // Graph for block 11
114  M_graph_block11.reset (new Epetra_FECrsGraph (Copy, * (M_pFESpace->map().map (Unique) ), 0) );
115  buildGraph ( elements (M_pFESpace->mesh() ),
116  quadRuleTetra4pt,
117  M_fespacePETA,
118  M_fespacePETA,
119  /*(14)*/ TAU_M*dot( grad(phi_i), grad(phi_j) )
120  ) >> M_graph_block11;
121  M_graph_block11->GlobalAssemble ( );
122  M_graph_block11->OptimizeStorage();
123 
124  }
125 
126  M_block_00.reset (new matrix_Type ( M_uFESpace->map(), *M_graph_block00 ) );
127  M_block_00->zero();
128  M_block_00->globalAssemble();
129 
130  M_block_01.reset (new matrix_Type ( M_uFESpace->map(), *M_graph_block01 ) );
131  M_block_01->zero();
132  M_block_01->globalAssemble( M_pFESpace->mapPtr(), M_uFESpace->mapPtr() );
133 
134  M_block_10.reset (new matrix_Type ( M_pFESpace->map(), *M_graph_block10 ) );
135  M_block_10->zero();
136  M_block_10->globalAssemble( M_uFESpace->mapPtr(), M_pFESpace->mapPtr() );
137 
138  M_block_11.reset (new matrix_Type ( M_pFESpace->map(), *M_graph_block11 ) );
139  M_block_11->zero();
140  M_block_11->globalAssemble( );
141 }
142 
143 void StabilizationSUPGALE::apply_matrix( const vector_Type& convective_velocity_previous_newton_step,
144  const vector_Type& velocity_previous_newton_step,
145  const vector_Type& pressure_previous_newton_step,
146  const vector_Type& velocity_rhs)
147 {
148  // missing force
149  if ( !M_useGraph )
150  {
151  M_block_00.reset (new matrix_Type ( M_uFESpace->map() ) );
152 
153  M_block_01.reset (new matrix_Type ( M_uFESpace->map() ) );
154 
155  M_block_10.reset (new matrix_Type ( M_pFESpace->map() ) );
156 
157  M_block_11.reset (new matrix_Type ( M_pFESpace->map() ) );
158  }
159 
160  M_block_00->zero();
161  M_block_01->zero();
162  M_block_10->zero();
163  M_block_11->zero();
164 
165  vector_Type beta_km1( convective_velocity_previous_newton_step, Repeated);
166  vector_Type u_km1( velocity_previous_newton_step, Repeated);
167  vector_Type p_km1( pressure_previous_newton_step, Repeated);
168  vector_Type u_bdf( velocity_rhs, Repeated);
169 
170  std::shared_ptr<SquareRoot_SUPGALE> squareroot(new SquareRoot_SUPGALE());
171 
172  MatrixSmall<3, 3> Eye;
173  Eye *= 0.0;
174  Eye[0][0] = 1;
175  Eye[1][1] = 1;
176  Eye[2][2] = 1;
177 
178  using namespace ExpressionAssembly;
179 
180  integrate(
181  elements(M_uFESpace->mesh()),
182  M_uFESpace->qr(),
183  M_fespaceUETA, // test w -> phi_i
184  M_fespaceUETA, // trial \delta u -> phi_j
185  /*(1)*/ TAU_M * value(M_density*M_density) * value(M_alpha/M_timestep) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), phi_j )
186  /*(2)*/ +TAU_M * value(M_density*M_density) * value(M_alpha/M_timestep) * dot( phi_j*grad(phi_i), value(M_fespaceUETA, u_km1) )
187  /*(3)*/ -TAU_M * value(M_density*M_density) * dot( phi_j*grad(phi_i), value(M_fespaceUETA, u_bdf) )
188  /*(4)*/ +TAU_M * value(M_density*M_density) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), phi_j*grad(M_fespaceUETA, u_km1) )
189  /*(5)*/ +TAU_M * value(M_density*M_density) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), value(M_fespaceUETA, beta_km1)*grad(phi_j) )
190  /*(6)*/ +TAU_M * value(M_density*M_density) * dot( phi_j*grad(phi_i), value(M_fespaceUETA, beta_km1)*grad(M_fespaceUETA, u_km1) )
191  /*(8)*/ +TAU_M * value(M_density) * dot( phi_j*grad(phi_i), grad(M_fespacePETA, p_km1) )
192  /*(9)*/ -TAU_M * value(M_density*M_viscosity) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), laplacian(phi_j) )
193  /*(10)*/ -TAU_M * value(M_density*M_viscosity) * dot( phi_j*grad(phi_i), laplacian(M_fespaceUETA, u_km1) )
194  /*(16)*/ +TAU_C * div(phi_i)*div(phi_j)
195  ) >> M_block_00;
196  M_block_00->globalAssemble();
197 
198  integrate(
199  elements(M_uFESpace->mesh()),
200  M_pFESpace->qr(),
201  M_fespacePETA, // test q -> phi_i
202  M_fespaceUETA, // trial \delta u -> phi_j
203  /*(11)*/ TAU_M * value(M_alpha*M_density/M_timestep) * dot( grad(phi_i), phi_j )
204  /*(12)*/ +TAU_M * value(M_density) * dot( grad(phi_i), grad(M_fespaceUETA, u_km1)*phi_j )
205  /*(13)*/ +TAU_M * value(M_density) * dot( grad(phi_i), value(M_fespaceUETA, beta_km1)*grad(phi_j) )
206  /*(15)*/ -TAU_M * value(M_viscosity) * dot( grad(phi_i), laplacian(phi_j) )
207  ) >> M_block_10;
208  M_block_10->globalAssemble( M_uFESpace->mapPtr(), M_pFESpace->mapPtr() );
209 
210  integrate(
211  elements(M_uFESpace->mesh()),
212  M_uFESpace->qr(),
213  M_fespaceUETA, // test w -> phi_i
214  M_fespacePETA, // trial \delta p -> phi_j
215  /*(7)*/ TAU_M * value(M_density) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), grad(phi_j) )
216  ) >> M_block_01;
217  M_block_01->globalAssemble( M_pFESpace->mapPtr(), M_uFESpace->mapPtr() );
218 
219  integrate(
220  elements(M_uFESpace->mesh()),
221  M_pFESpace->qr(),
222  M_fespacePETA, // test q -> phi_i
223  M_fespacePETA, // trial \delta p -> phi_j
224  /*(14)*/ TAU_M*dot( grad(phi_i), grad(phi_j) )
225  ) >> M_block_11;
226  M_block_11->globalAssemble();
227 
228 }
229 
230 void StabilizationSUPGALE::apply_vector( vectorPtr_Type& residual_velocity,
231  vectorPtr_Type& residual_pressure,
232  const vector_Type& convective_velocity_previous_newton_step,
233  const vector_Type& velocity_previous_newton_step,
234  const vector_Type& pressure_previous_newton_step,
235  const vector_Type& velocity_rhs)
236 {
237  // missing force, terms 11 and 12
238 
239  vector_Type beta_km1( convective_velocity_previous_newton_step, Repeated);
240  vector_Type u_km1( velocity_previous_newton_step, Repeated);
241  vector_Type p_km1( pressure_previous_newton_step, Repeated);
242  vector_Type u_bdf( velocity_rhs, Repeated);
243 
244  std::shared_ptr<SquareRoot_SUPGALE> squareroot(new SquareRoot_SUPGALE());
245 
246  // Matrix needed to evaluate the divergence of a vector (term with TAU_C)
247  MatrixSmall<3, 3> Eye;
248  Eye *= 0.0;
249  Eye[0][0] = 1;
250  Eye[1][1] = 1;
251  Eye[2][2] = 1;
252 
253  using namespace ExpressionAssembly;
254 
255  integrate(
256  elements(M_uFESpace->mesh()),
257  M_uFESpace->qr(),
258  M_fespaceUETA,
259  /*(1)*/ TAU_M * value(M_density*M_density)*value(M_alpha/M_timestep) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), value(M_fespaceUETA, u_km1) )
260  /*(2)*/ -TAU_M * value(M_density*M_density) * dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), value(M_fespaceUETA, u_bdf) )
261  /*(3)*/ +TAU_M*value(M_density*M_density)*dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), value(M_fespaceUETA, beta_km1)*grad(M_fespaceUETA, u_km1) )
262  /*(4)*/ +TAU_M*value(M_density)*dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), grad(M_fespacePETA, p_km1) )
263  /*(5)*/ -TAU_M*value(M_density*M_viscosity)*dot( value(M_fespaceUETA, beta_km1)*grad(phi_i), laplacian(M_fespaceUETA, u_km1) )
264  /*(11)*/ +TAU_C*div(phi_i)* trace ( grad(M_fespaceUETA, u_km1) )
265  )
266  >> residual_velocity;
267 
268  integrate(
269  elements(M_uFESpace->mesh()),
270  M_pFESpace->qr(),
271  M_fespacePETA,
272  /*(6)*/ TAU_M*value(M_density*M_alpha/M_timestep)*dot( grad(phi_i), value(M_fespaceUETA, u_km1) )
273  /*(7)*/ -TAU_M*value(M_density)*dot( grad(phi_i), value(M_fespaceUETA, u_bdf))
274  /*(8)*/ +TAU_M*value(M_density)*dot( grad(phi_i), value(M_fespaceUETA, beta_km1)*grad(M_fespaceUETA, u_km1) )
275  /*(9)*/ +TAU_M*dot( grad(phi_i), grad(M_fespacePETA,p_km1) )
276  /*(10)*/ -TAU_M*value(M_viscosity)*dot(grad(phi_i), laplacian(M_fespaceUETA, u_km1))
277  )
278  >> residual_pressure;
279 }
280 
281 } // namespace LifeV
TAU_M_DEN
#define TAU_M_DEN
Definition: StabilizationSUPG.cpp:5
TAU_M_DEN_VEL
#define TAU_M_DEN_VEL
Definition: StabilizationSUPG.cpp:7
ASSERT
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90
TAU_M_DEN_VISC
#define TAU_M_DEN_VISC
Definition: StabilizationSUPG.cpp:8
TAU_C
#define TAU_C
Definition: StabilizationSUPG.cpp:11
TAU_M
#define TAU_M
Definition: StabilizationSUPG.cpp:4
TAU_M_DEN_DT
#define TAU_M_DEN_DT
Definition: StabilizationSUPG.cpp:6