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