ExaDG
Loading...
Searching...
No Matches
weak_boundary_conditions.h
1/*
2 * weak_boundary_conditions.h
3 *
4 * Created on: Jun 13, 2019
5 * Author: fehn
6 */
7
8#ifndef INCLUDE_EXADG_INCOMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_OPERATORS_WEAK_BOUNDARY_CONDITIONS_H_
9#define INCLUDE_EXADG_INCOMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_OPERATORS_WEAK_BOUNDARY_CONDITIONS_H_
10
11#include <exadg/functions_and_boundary_conditions/evaluate_functions.h>
12#include <exadg/incompressible_navier_stokes/user_interface/boundary_descriptor.h>
13#include <exadg/incompressible_navier_stokes/user_interface/enum_types.h>
14#include <exadg/matrix_free/integrators.h>
15#include <exadg/operators/operator_type.h>
16
17namespace ExaDG
18{
19namespace IncNS
20{
21// clang-format off
22/*
23 * Velocity:
24 *
25 * The following two functions calculate the interior/exterior value for boundary faces depending on the
26 * operator type, the type of the boundary face and the given boundary conditions.
27 *
28 * +-------------------------+--------------------+------------------------------+
29 * | Dirichlet boundaries | Neumann boundaries | symmetry boundaries |
30 * +-------------------------+-------------------------+--------------------+------------------------------+
31 * | full operator | u⁺ = -u⁻ + 2g | u⁺ = u⁻ | u⁺ = u⁻ - 2 (u⁻*n)n |
32 * +-------------------------+-------------------------+--------------------+------------------------------+
33 * | homogeneous operator | u⁺ = -u⁻ | u⁺ = u⁻ | u⁺ = u⁻ - 2 (u⁻*n)n |
34 * +-------------------------+-------------------------+--------------------+------------------------------+
35 * | inhomogeneous operator | u⁺ = -u⁻ + 2g , u⁻ = 0 | u⁺ = u⁻ , u⁻ = 0 | u⁺ = u⁻ - 2 (u⁻*n)n , u⁻ = 0 |
36 * +-------------------------+-------------------------+--------------------+------------------------------+
37 *
38 */
39// clang-format on
40template<int dim, typename Number>
41inline DEAL_II_ALWAYS_INLINE //
42 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
43 calculate_interior_value(unsigned int const q,
44 FaceIntegrator<dim, dim, Number> const & integrator,
45 OperatorType const & operator_type)
46{
47 // element e⁻
48 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> value_m;
49
50 if(operator_type == OperatorType::full or operator_type == OperatorType::homogeneous)
51 {
52 value_m = integrator.get_value(q);
53 }
54 else if(operator_type == OperatorType::inhomogeneous)
55 {
56 // do nothing, value_m is already initialized with zeros
57 }
58 else
59 {
60 AssertThrow(false, dealii::ExcMessage("Specified OperatorType is not implemented!"));
61 }
62
63 return value_m;
64}
65
66template<int dim, typename Number>
67inline DEAL_II_ALWAYS_INLINE //
68 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
69 calculate_exterior_value(dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & value_m,
70 unsigned int const q,
71 FaceIntegrator<dim, dim, Number> const & integrator,
72 OperatorType const & operator_type,
73 BoundaryTypeU const & boundary_type,
74 dealii::types::boundary_id const boundary_id,
75 std::shared_ptr<BoundaryDescriptorU<dim> const> boundary_descriptor,
76 double const & time)
77{
78 // element e⁺
79 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> value_p;
80
81 if(boundary_type == BoundaryTypeU::Dirichlet or boundary_type == BoundaryTypeU::DirichletCached)
82 {
83 if(operator_type == OperatorType::full or operator_type == OperatorType::inhomogeneous)
84 {
85 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> g;
86
87 if(boundary_type == BoundaryTypeU::Dirichlet)
88 {
89 auto bc = boundary_descriptor->dirichlet_bc.find(boundary_id)->second;
90 auto q_points = integrator.quadrature_point(q);
91
92 g = FunctionEvaluator<1, dim, Number>::value(*bc, q_points, time);
93 }
94 else if(boundary_type == BoundaryTypeU::DirichletCached)
95 {
96 auto bc = boundary_descriptor->get_dirichlet_cached_data();
97 g = FunctionEvaluator<1, dim, Number>::value(*bc,
98 integrator.get_current_cell_index(),
99 q,
100 integrator.get_quadrature_index());
101 }
102 else
103 {
104 AssertThrow(false, dealii::ExcMessage("Not implemented."));
105 }
106
107 value_p = -value_m + dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>(2.0 * g);
108 }
109 else if(operator_type == OperatorType::homogeneous)
110 {
111 value_p = -value_m;
112 }
113 else
114 {
115 AssertThrow(false, dealii::ExcMessage("Specified OperatorType is not implemented!"));
116 }
117 }
118 else if(boundary_type == BoundaryTypeU::Neumann)
119 {
120 value_p = value_m;
121 }
122 else if(boundary_type == BoundaryTypeU::Symmetry)
123 {
124 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> normal_m =
125 integrator.get_normal_vector(q);
126
127 value_p = value_m - 2.0 * (value_m * normal_m) * normal_m;
128 }
129 else
130 {
131 AssertThrow(false, dealii::ExcMessage("Boundary type of face is invalid or not implemented."));
132 }
133
134 return value_p;
135}
136
137/*
138 * This function calculates the exterior velocity on boundary faces
139 * according to:
140 *
141 * Dirichlet boundary: u⁺ = -u⁻ + 2g (mirror) or g (direct)
142 * Neumann boundary: u⁺ = u⁻
143 * symmetry boundary: u⁺ = u⁻ -(u⁻*n)n - (u⁻*n)n = u⁻ - 2 (u⁻*n)n
144 *
145 * The name "nonlinear" indicates that this function is used when
146 * evaluating the nonlinear convective operator.
147 */
148template<int dim, typename Number>
149inline DEAL_II_ALWAYS_INLINE //
150 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
151 calculate_exterior_value_nonlinear(
152 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & u_m,
153 unsigned int const q,
154 FaceIntegrator<dim, dim, Number> & integrator,
155 BoundaryTypeU const & boundary_type,
156 TypeDirichletBCs const & type_dirichlet_bc,
157 dealii::types::boundary_id const boundary_id,
158 std::shared_ptr<BoundaryDescriptorU<dim> const> boundary_descriptor,
159 double const & time)
160{
161 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> u_p;
162
163 if(boundary_type == BoundaryTypeU::Dirichlet or boundary_type == BoundaryTypeU::DirichletCached)
164 {
165 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> g;
166
167 if(boundary_type == BoundaryTypeU::Dirichlet)
168 {
169 auto bc = boundary_descriptor->dirichlet_bc.find(boundary_id)->second;
170 auto q_points = integrator.quadrature_point(q);
171
172 g = FunctionEvaluator<1, dim, Number>::value(*bc, q_points, time);
173 }
174 else if(boundary_type == BoundaryTypeU::DirichletCached)
175 {
176 auto bc = boundary_descriptor->get_dirichlet_cached_data();
177 g = FunctionEvaluator<1, dim, Number>::value(*bc,
178 integrator.get_current_cell_index(),
179 q,
180 integrator.get_quadrature_index());
181 }
182 else
183 {
184 AssertThrow(false, dealii::ExcMessage("Not implemented."));
185 }
186
187 if(type_dirichlet_bc == TypeDirichletBCs::Mirror)
188 {
189 u_p = -u_m + dealii::make_vectorized_array<Number>(2.0) * g;
190 }
191 else if(type_dirichlet_bc == TypeDirichletBCs::Direct)
192 {
193 u_p = g;
194 }
195 else
196 {
197 AssertThrow(false, dealii::ExcMessage("TypeDirichletBCs is not implemented."));
198 }
199 }
200 else if(boundary_type == BoundaryTypeU::Neumann)
201 {
202 u_p = u_m;
203 }
204 else if(boundary_type == BoundaryTypeU::Symmetry)
205 {
206 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> normal_m =
207 integrator.get_normal_vector(q);
208
209 u_p = u_m - 2. * (u_m * normal_m) * normal_m;
210 }
211 else
212 {
213 AssertThrow(false, dealii::ExcMessage("Boundary type of face is invalid or not implemented."));
214 }
215
216 return u_p;
217}
218
219template<int dim, typename Number>
220inline DEAL_II_ALWAYS_INLINE //
221 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
222 calculate_exterior_value_from_dof_vector(
223 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & value_m,
224 unsigned int const q,
225 FaceIntegrator<dim, dim, Number> const & integrator_bc,
226 OperatorType const & operator_type,
227 BoundaryTypeU const & boundary_type)
228{
229 // element e⁺
230 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> value_p;
231
232 if(boundary_type == BoundaryTypeU::Dirichlet or boundary_type == BoundaryTypeU::DirichletCached)
233 {
234 if(operator_type == OperatorType::full or operator_type == OperatorType::inhomogeneous)
235 {
236 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> g = integrator_bc.get_value(q);
237
238 value_p = -value_m + dealii::make_vectorized_array<Number>(2.0) * g;
239 }
240 else if(operator_type == OperatorType::homogeneous)
241 {
242 value_p = -value_m;
243 }
244 else
245 {
246 AssertThrow(false, dealii::ExcMessage("Specified OperatorType is not implemented!"));
247 }
248 }
249 else if(boundary_type == BoundaryTypeU::Neumann)
250 {
251 value_p = value_m;
252 }
253 else if(boundary_type == BoundaryTypeU::Symmetry)
254 {
255 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> normal_m =
256 integrator_bc.get_normal_vector(q);
257
258 value_p = value_m - 2.0 * (value_m * normal_m) * normal_m;
259 }
260 else
261 {
262 AssertThrow(false, dealii::ExcMessage("Boundary type of face is invalid or not implemented."));
263 }
264
265 return value_p;
266}
267
268/*
269 * Pressure:
270 *
271 * These two function calculate the interior/exterior value on boundary faces depending on the
272 * operator type, the type of the boundary face and the given boundary conditions.
273 *
274 * +--------------------+----------------------+
275 * | Neumann boundaries | Dirichlet boundaries |
276 * +-------------------------+--------------------+----------------------+
277 * | full operator | p⁺ = p⁻ | p⁺ = - p⁻ + 2g |
278 * +-------------------------+--------------------+----------------------+
279 * | homogeneous operator | p⁺ = p⁻ | p⁺ = - p⁻ |
280 * +-------------------------+--------------------+----------------------+
281 * | inhomogeneous operator | p⁺ = 0 , p⁻ = 0 | p⁺ = 2g , p⁻ = 0 |
282 * +-------------------------+--------------------+----------------------+
283 *
284 */
285template<int dim, typename Number>
286inline DEAL_II_ALWAYS_INLINE //
287 dealii::VectorizedArray<Number>
288 calculate_interior_value(unsigned int const q,
289 FaceIntegrator<dim, 1, Number> const & integrator,
290 OperatorType const & operator_type)
291{
292 // element e⁻
293 dealii::VectorizedArray<Number> value_m = dealii::make_vectorized_array<Number>(0.0);
294
295 if(operator_type == OperatorType::full or operator_type == OperatorType::homogeneous)
296 {
297 value_m = integrator.get_value(q);
298 }
299 else if(operator_type == OperatorType::inhomogeneous)
300 {
301 // do nothing, value_m is already initialized with zeros
302 }
303 else
304 {
305 AssertThrow(false, dealii::ExcMessage("Specified OperatorType is not implemented!"));
306 }
307
308 return value_m;
309}
310
311
312template<int dim, typename Number>
313inline DEAL_II_ALWAYS_INLINE //
314 dealii::VectorizedArray<Number>
315 calculate_exterior_value(dealii::VectorizedArray<Number> const & value_m,
316 unsigned int const q,
317 FaceIntegrator<dim, 1, Number> const & integrator,
318 OperatorType const & operator_type,
319 BoundaryTypeP const & boundary_type,
320 dealii::types::boundary_id const boundary_id,
321 std::shared_ptr<BoundaryDescriptorP<dim> const> boundary_descriptor,
322 double const & time,
323 double const & inverse_scaling_factor)
324{
325 dealii::VectorizedArray<Number> value_p = dealii::make_vectorized_array<Number>(0.0);
326
327 if(boundary_type == BoundaryTypeP::Dirichlet)
328 {
329 if(operator_type == OperatorType::full or operator_type == OperatorType::inhomogeneous)
330 {
331 auto bc = boundary_descriptor->dirichlet_bc.find(boundary_id)->second;
332 auto q_points = integrator.quadrature_point(q);
333
334 dealii::VectorizedArray<Number> g =
335 FunctionEvaluator<0, dim, Number>::value(*bc, q_points, time);
336
337 value_p = -value_m + 2.0 * inverse_scaling_factor * g;
338 }
339 else if(operator_type == OperatorType::homogeneous)
340 {
341 value_p = -value_m;
342 }
343 else
344 {
345 AssertThrow(false, dealii::ExcMessage("Specified OperatorType is not implemented!"));
346 }
347 }
348 else if(boundary_type == BoundaryTypeP::Neumann)
349 {
350 value_p = value_m;
351 }
352 else
353 {
354 AssertThrow(false, dealii::ExcMessage("Boundary type of face is invalid or not implemented."));
355 }
356
357 return value_p;
358}
359
360template<int dim, typename Number>
361inline DEAL_II_ALWAYS_INLINE //
362 dealii::VectorizedArray<Number>
363 calculate_exterior_value_from_dof_vector(dealii::VectorizedArray<Number> const & value_m,
364 unsigned int const q,
365 FaceIntegrator<dim, 1, Number> const & integrator_bc,
366 OperatorType const & operator_type,
367 BoundaryTypeP const & boundary_type,
368 double const & inverse_scaling_factor)
369{
370 dealii::VectorizedArray<Number> value_p = dealii::make_vectorized_array<Number>(0.0);
371
372 if(boundary_type == BoundaryTypeP::Dirichlet)
373 {
374 if(operator_type == OperatorType::full or operator_type == OperatorType::inhomogeneous)
375 {
376 dealii::VectorizedArray<Number> g = integrator_bc.get_value(q);
377
378 value_p = -value_m + 2.0 * inverse_scaling_factor * g;
379 }
380 else if(operator_type == OperatorType::homogeneous)
381 {
382 value_p = -value_m;
383 }
384 else
385 {
386 AssertThrow(false, dealii::ExcMessage("Specified OperatorType is not implemented!"));
387 }
388 }
389 else if(boundary_type == BoundaryTypeP::Neumann)
390 {
391 value_p = value_m;
392 }
393 else
394 {
395 AssertThrow(false, dealii::ExcMessage("Boundary type of face is invalid or not implemented."));
396 }
397
398 return value_p;
399}
400
401// clang-format off
402/*
403 * This function calculates the exterior velocity gradient in normal direction depending
404 * on the operator type, the type of the boundary face and the given boundary conditions.
405 *
406 * The function argument normal_gradient_m has the meaning of F(u⁻)*n with
407 *
408 * Divergence formulation: F(u) = F_nu(u) / nu = ( grad(u) + grad(u)^T )
409 * Laplace formulation: F(u) = F_nu(u) / nu = grad(u)
410 *
411 * +---------------------------------+---------------------------------------+----------------------------------------------------+
412 * | Dirichlet boundaries | Neumann boundaries | symmetry boundaries |
413 * +-------------------------+---------------------------------+---------------------------------------+----------------------------------------------------+
414 * | full operator | F(u⁺)*n = F(u⁻)*n | F(u⁺)*n = -F(u⁻)*n + 2h | F(u⁺)*n = -F(u⁻)*n + 2*[(F(u⁻)*n)*n]n |
415 * +-------------------------+---------------------------------+---------------------------------------+----------------------------------------------------+
416 * | homogeneous operator | F(u⁺)*n = F(u⁻)*n | F(u⁺)*n = -F(u⁻)*n | F(u⁺)*n = -F(u⁻)*n + 2*[(F(u⁻)*n)*n]n |
417 * +-------------------------+---------------------------------+---------------------------------------+----------------------------------------------------+
418 * | inhomogeneous operator | F(u⁺)*n = F(u⁻)*n, F(u⁻)*n = 0 | F(u⁺)*n = -F(u⁻)*n + 2h , F(u⁻)*n = 0 | F(u⁺)*n = -F(u⁻)*n + 2*[(F(u⁻)*n)*n]n, F(u⁻)*n = 0 |
419 * +-------------------------+---------------------------------+---------------------------------------+----------------------------------------------------+
420 *
421 * +---------------------------------+---------------------------------------+----------------------------------------------------+
422 * | Dirichlet boundaries | Neumann boundaries | symmetry boundaries |
423 * +-------------------------+---------------------------------+---------------------------------------+----------------------------------------------------+
424 * | full operator | {{F(u)}}*n = F(u⁻)*n | {{F(u)}}*n = h | {{F(u)}}*n = [(F(u⁻)*n)*n]n |
425 * +-------------------------+---------------------------------+---------------------------------------+----------------------------------------------------+
426 * | homogeneous operator | {{F(u)}}*n = F(u⁻)*n | {{F(u)}}*n = 0 | {{F(u)}}*n = [(F(u⁻)*n)*n]n |
427 * +-------------------------+---------------------------------+---------------------------------------+----------------------------------------------------+
428 * | inhomogeneous operator | {{F(u)}}*n = 0 | {{F(u)}}*n = h | {{F(u)}}*n = 0 |
429 * +-------------------------+---------------------------------+---------------------------------------+----------------------------------------------------+
430 */
431// clang-format on
432
433template<int dim, typename Number>
434inline DEAL_II_ALWAYS_INLINE //
435 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
436 calculate_exterior_normal_gradient(
437 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> const & normal_gradient_m,
438 unsigned int const q,
439 FaceIntegrator<dim, dim, Number> const & integrator,
440 OperatorType const & operator_type,
441 BoundaryTypeU const & boundary_type,
442 dealii::types::boundary_id const boundary_id,
443 std::shared_ptr<BoundaryDescriptorU<dim> const> boundary_descriptor,
444 double const & time,
445 bool const variable_normal_vector)
446{
447 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> normal_gradient_p;
448
449 if(boundary_type == BoundaryTypeU::Dirichlet or boundary_type == BoundaryTypeU::DirichletCached)
450 {
451 normal_gradient_p = normal_gradient_m;
452 }
453 else if(boundary_type == BoundaryTypeU::Neumann)
454 {
455 if(operator_type == OperatorType::full or operator_type == OperatorType::inhomogeneous)
456 {
457 auto bc = boundary_descriptor->neumann_bc.find(boundary_id)->second;
458 auto q_points = integrator.quadrature_point(q);
459
460 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> h;
461 if(variable_normal_vector == false)
462 {
463 h = FunctionEvaluator<1, dim, Number>::value(*bc, q_points, time);
464 }
465 else
466 {
467 auto normals_m = integrator.get_normal_vector(q);
468 h = FunctionEvaluator<1, dim, Number>::value(
469 *(std::dynamic_pointer_cast<FunctionWithNormal<dim>>(bc)), q_points, normals_m, time);
470 }
471
472 normal_gradient_p = -normal_gradient_m + 2.0 * h;
473 }
474 else if(operator_type == OperatorType::homogeneous)
475 {
476 normal_gradient_p = -normal_gradient_m;
477 }
478 else
479 {
480 AssertThrow(false, dealii::ExcMessage("Specified OperatorType is not implemented!"));
481 }
482 }
483 else if(boundary_type == BoundaryTypeU::Symmetry)
484 {
485 auto normal_m = integrator.get_normal_vector(q);
486 normal_gradient_p = -normal_gradient_m + 2.0 * (normal_gradient_m * normal_m) * normal_m;
487 }
488 else
489 {
490 AssertThrow(false, dealii::ExcMessage("Boundary type of face is invalid or not implemented."));
491 }
492
493 return normal_gradient_p;
494}
495
496} // namespace IncNS
497} // namespace ExaDG
498
499#endif /* INCLUDE_EXADG_INCOMPRESSIBLE_NAVIER_STOKES_SPATIAL_DISCRETIZATION_OPERATORS_WEAK_BOUNDARY_CONDITIONS_H_ \
500 */
ExaDG
Definition driver.cpp:33
Generated by doxygen 1.10.0