52 typedef dealii::VectorizedArray<Number> scalar;
53 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
54 typedef dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> tensor;
56 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
58 typedef CellIntegrator<dim, dim, Number> IntegratorCell;
59 typedef FaceIntegrator<dim, dim, Number> IntegratorFace;
63 reinit(dealii::MatrixFree<dim, Number>
const & matrix_free,
65 unsigned int const dof_index,
66 unsigned int const quad_index_linearized,
67 bool const use_own_velocity_storage)
73 std::make_shared<IntegratorCell>(matrix_free, dof_index, quad_index_linearized);
74 integrator_velocity_m =
75 std::make_shared<IntegratorFace>(matrix_free,
true, dof_index, quad_index_linearized);
76 integrator_velocity_p =
77 std::make_shared<IntegratorFace>(matrix_free,
false, dof_index, quad_index_linearized);
81 integrator_grid_velocity =
82 std::make_shared<IntegratorCell>(matrix_free, dof_index, quad_index_linearized);
83 integrator_grid_velocity_face =
84 std::make_shared<IntegratorFace>(matrix_free,
true, dof_index, quad_index_linearized);
87 if(use_own_velocity_storage)
90 matrix_free.initialize_dof_vector(velocity.own(), dof_index);
95 matrix_free.initialize_dof_vector(grid_velocity.own(), dof_index);
97 AssertThrow(data.formulation == FormulationConvectiveTerm::ConvectiveFormulation,
99 "ALE formulation can only be used in combination with ConvectiveFormulation"));
108 flags.cells = dealii::update_JxW_values | dealii::update_gradients;
109 flags.inner_faces = dealii::update_JxW_values | dealii::update_normal_vectors;
110 flags.boundary_faces = dealii::update_JxW_values | dealii::update_normal_vectors;
120 get_integrator_flags()
const
124 if(data.formulation == FormulationConvectiveTerm::DivergenceFormulation)
126 flags.cell_evaluate = dealii::EvaluationFlags::values;
127 flags.cell_integrate = dealii::EvaluationFlags::gradients;
129 flags.face_evaluate = dealii::EvaluationFlags::values;
130 flags.face_integrate = dealii::EvaluationFlags::values;
132 else if(data.formulation == FormulationConvectiveTerm::ConvectiveFormulation)
134 flags.cell_evaluate = dealii::EvaluationFlags::values | dealii::EvaluationFlags::gradients;
135 flags.cell_integrate = dealii::EvaluationFlags::values;
137 flags.face_evaluate = dealii::EvaluationFlags::values;
138 flags.face_integrate = dealii::EvaluationFlags::values;
142 AssertThrow(
false, dealii::ExcMessage(
"Not implemented."));
161 set_velocity_copy(VectorType
const & src)
163 velocity.own() = src;
165 velocity->update_ghost_values();
169 set_velocity_ptr(VectorType
const & src)
173 velocity->update_ghost_values();
177 set_grid_velocity_ptr(VectorType
const & src)
179 grid_velocity.reset(src);
181 grid_velocity->update_ghost_values();
185 get_grid_velocity()
const
187 return *grid_velocity;
190 inline DEAL_II_ALWAYS_INLINE
192 get_velocity_cell(
unsigned int const q)
const
194 return integrator_velocity->get_value(q);
197 inline DEAL_II_ALWAYS_INLINE
199 get_velocity_gradient_cell(
unsigned int const q)
const
201 return integrator_velocity->get_gradient(q);
204 inline DEAL_II_ALWAYS_INLINE
206 get_velocity_m(
unsigned int const q)
const
208 return integrator_velocity_m->get_value(q);
211 inline DEAL_II_ALWAYS_INLINE
213 get_velocity_p(
unsigned int const q)
const
215 return integrator_velocity_p->get_value(q);
219 inline DEAL_II_ALWAYS_INLINE
221 get_grid_velocity_cell(
unsigned int const q)
const
223 return integrator_grid_velocity->get_value(q);
229 inline DEAL_II_ALWAYS_INLINE
231 get_grid_velocity_face(
unsigned int const q)
const
233 return integrator_grid_velocity_face->get_value(q);
238 reinit_cell(
unsigned int const cell)
const
240 integrator_velocity->reinit(cell);
243 integrator_grid_velocity->reinit(cell);
245 if(data.formulation == FormulationConvectiveTerm::DivergenceFormulation)
247 integrator_velocity->gather_evaluate(*velocity, dealii::EvaluationFlags::values);
249 else if(data.formulation == FormulationConvectiveTerm::ConvectiveFormulation)
251 integrator_velocity->gather_evaluate(*velocity,
252 dealii::EvaluationFlags::values |
253 dealii::EvaluationFlags::gradients);
256 integrator_grid_velocity->gather_evaluate(*grid_velocity, dealii::EvaluationFlags::values);
260 AssertThrow(
false, dealii::ExcMessage(
"Not implemented."));
265 reinit_face(
unsigned int const face)
const
267 integrator_velocity_m->reinit(face);
268 integrator_velocity_m->gather_evaluate(*velocity, dealii::EvaluationFlags::values);
270 integrator_velocity_p->reinit(face);
271 integrator_velocity_p->gather_evaluate(*velocity, dealii::EvaluationFlags::values);
275 integrator_grid_velocity_face->reinit(face);
276 integrator_grid_velocity_face->gather_evaluate(*grid_velocity,
277 dealii::EvaluationFlags::values);
282 reinit_boundary_face(
unsigned int const face)
const
284 integrator_velocity_m->reinit(face);
285 integrator_velocity_m->gather_evaluate(*velocity, dealii::EvaluationFlags::values);
289 integrator_grid_velocity_face->reinit(face);
290 integrator_grid_velocity_face->gather_evaluate(*grid_velocity,
291 dealii::EvaluationFlags::values);
296 reinit_face_cell_based(
unsigned int const cell,
297 unsigned int const face,
298 dealii::types::boundary_id
const boundary_id)
const
300 integrator_velocity_m->reinit(cell, face);
301 integrator_velocity_m->gather_evaluate(*velocity, dealii::EvaluationFlags::values);
305 integrator_grid_velocity_face->reinit(cell, face);
306 integrator_grid_velocity_face->gather_evaluate(*grid_velocity,
307 dealii::EvaluationFlags::values);
310 if(boundary_id == dealii::numbers::internal_face_boundary_id)
324 inline DEAL_II_ALWAYS_INLINE
326 get_volume_flux_linearized_divergence_formulation(vector
const & delta_u,
327 unsigned int const q)
const
329 vector u = get_velocity_cell(q);
330 tensor F = outer_product(u, delta_u);
333 return -(F + transpose(F));
336 inline DEAL_II_ALWAYS_INLINE
338 get_volume_flux_linearized_convective_formulation(vector
const & delta_u,
339 tensor
const & grad_delta_u,
340 unsigned int const q)
const
343 vector w = get_velocity_cell(q);
344 tensor grad_u = get_velocity_gradient_cell(q);
348 w -= get_grid_velocity_cell(q);
352 vector F = grad_u * delta_u + grad_delta_u * w;
363 inline DEAL_II_ALWAYS_INLINE
364 std::tuple<vector, vector>
365 calculate_flux_nonlinear_interior_and_neighbor(vector
const & uM,
367 vector
const & normalM,
368 vector
const & u_grid)
const
370 vector flux_m, flux_p;
372 if(data.formulation == FormulationConvectiveTerm::DivergenceFormulation)
374 vector flux = calculate_lax_friedrichs_flux(uM, uP, normalM);
379 else if(data.formulation == FormulationConvectiveTerm::ConvectiveFormulation)
382 scalar average_u_normal = 0.5 * (uM + uP) * normalM;
384 average_u_normal -= u_grid * normalM;
386 flux = calculate_upwind_flux(uM, uP, average_u_normal);
390 flux_m = flux - average_u_normal * uM;
391 flux_p = -flux + average_u_normal * uP;
395 AssertThrow(
false, dealii::ExcMessage(
"Not implemented."));
398 return std::make_tuple(flux_m, flux_p);
406 inline DEAL_II_ALWAYS_INLINE
408 calculate_flux_nonlinear_boundary(vector
const & uM,
410 vector
const & normalM,
411 vector
const & u_grid,
412 BoundaryTypeU
const & boundary_type)
const
416 if(data.formulation == FormulationConvectiveTerm::DivergenceFormulation)
418 flux = calculate_lax_friedrichs_flux(uM, uP, normalM);
420 if(boundary_type == BoundaryTypeU::Neumann and data.use_outflow_bc ==
true)
421 apply_outflow_bc(flux, uM * normalM);
423 else if(data.formulation == FormulationConvectiveTerm::ConvectiveFormulation)
425 scalar average_u_normal = 0.5 * (uM + uP) * normalM;
427 average_u_normal -= u_grid * normalM;
429 flux = calculate_upwind_flux(uM, uP, average_u_normal);
431 if(boundary_type == BoundaryTypeU::Neumann and data.use_outflow_bc ==
true)
433 if(data.ale ==
false)
434 apply_outflow_bc(flux, uM * normalM);
436 apply_outflow_bc(flux, (uM - u_grid) * normalM);
441 flux = flux - average_u_normal * uM;
445 AssertThrow(
false, dealii::ExcMessage(
"Not implemented."));
457 inline DEAL_II_ALWAYS_INLINE
458 std::tuple<vector, vector>
459 calculate_flux_linearized_interior_and_neighbor(vector
const & uM,
461 vector
const & delta_uM,
462 vector
const & delta_uP,
463 vector
const & normalM,
464 unsigned int const q)
const
468 if(data.formulation == FormulationConvectiveTerm::DivergenceFormulation)
470 fluxM = calculate_lax_friedrichs_flux_linearized(uM, uP, delta_uM, delta_uP, normalM);
473 else if(data.formulation == FormulationConvectiveTerm::ConvectiveFormulation)
477 u_grid = get_grid_velocity_face(q);
479 vector flux = calculate_upwind_flux_linearized(uM, uP, u_grid, delta_uM, delta_uP, normalM);
481 scalar average_u_normal = 0.5 * (uM + uP) * normalM;
483 average_u_normal -= u_grid * normalM;
485 scalar average_delta_u_normal = 0.5 * (delta_uM + delta_uP) * normalM;
489 fluxM = flux - average_u_normal * delta_uM - average_delta_u_normal * uM;
491 fluxP = -flux + average_u_normal * delta_uP + average_delta_u_normal * uP;
495 AssertThrow(
false, dealii::ExcMessage(
"Not implemented."));
498 return std::make_tuple(fluxM, fluxP);
505 inline DEAL_II_ALWAYS_INLINE
507 calculate_flux_linearized_interior(vector
const & uM,
509 vector
const & delta_uM,
510 vector
const & delta_uP,
511 vector
const & normalM,
512 unsigned int const q)
const
516 if(data.formulation == FormulationConvectiveTerm::DivergenceFormulation)
518 flux = calculate_lax_friedrichs_flux_linearized(uM, uP, delta_uM, delta_uP, normalM);
520 else if(data.formulation == FormulationConvectiveTerm::ConvectiveFormulation)
524 u_grid = get_grid_velocity_face(q);
526 flux = calculate_upwind_flux_linearized(uM, uP, u_grid, delta_uM, delta_uP, normalM);
528 scalar average_u_normal = 0.5 * (uM + uP) * normalM;
531 average_u_normal -= u_grid * normalM;
533 scalar average_delta_u_normal = 0.5 * (delta_uM + delta_uP) * normalM;
537 flux = flux - average_u_normal * delta_uM - average_delta_u_normal * uM;
541 AssertThrow(
false, dealii::ExcMessage(
"Not implemented."));
553 inline DEAL_II_ALWAYS_INLINE
555 calculate_flux_linearized_boundary(vector
const & uM,
557 vector
const & delta_uM,
558 vector
const & delta_uP,
559 vector
const & normalM,
560 BoundaryTypeU
const & boundary_type,
561 unsigned int const q)
const
565 if(data.formulation == FormulationConvectiveTerm::DivergenceFormulation)
567 flux = calculate_lax_friedrichs_flux_linearized(uM, uP, delta_uM, delta_uP, normalM);
569 if(boundary_type == BoundaryTypeU::Neumann and data.use_outflow_bc ==
true)
570 apply_outflow_bc(flux, uM * normalM);
572 else if(data.formulation == FormulationConvectiveTerm::ConvectiveFormulation)
576 u_grid = get_grid_velocity_face(q);
578 flux = calculate_upwind_flux_linearized(uM, uP, u_grid, delta_uM, delta_uP, normalM);
580 if(boundary_type == BoundaryTypeU::Neumann and data.use_outflow_bc ==
true)
581 apply_outflow_bc(flux, uM * normalM);
583 scalar average_u_normal = 0.5 * (uM + uP) * normalM;
585 average_u_normal -= u_grid * normalM;
587 scalar average_delta_u_normal = 0.5 * (delta_uM + delta_uP) * normalM;
591 flux = flux - average_u_normal * delta_uM - average_delta_u_normal * uM;
595 AssertThrow(
false, dealii::ExcMessage(
"Not implemented."));
607 inline DEAL_II_ALWAYS_INLINE
609 calculate_lambda(scalar
const & uM_n, scalar
const & uP_n)
const
611 return data.upwind_factor * 2.0 * std::max(std::abs(uM_n), std::abs(uP_n));
617 inline DEAL_II_ALWAYS_INLINE
619 calculate_lax_friedrichs_flux(vector
const & uM,
621 vector
const & normalM)
const
623 scalar uM_n = uM * normalM;
624 scalar uP_n = uP * normalM;
626 vector average_normal_flux =
627 dealii::make_vectorized_array<Number>(0.5) * (uM * uM_n + uP * uP_n);
629 vector jump_value = uM - uP;
631 scalar lambda = calculate_lambda(uM_n, uP_n);
633 return (average_normal_flux + 0.5 * lambda * jump_value);
639 inline DEAL_II_ALWAYS_INLINE
641 calculate_lax_friedrichs_flux_linear_transport(vector
const & uM,
645 vector
const & normalM)
const
647 scalar wM_n = wM * normalM;
648 scalar wP_n = wP * normalM;
650 vector average_normal_flux =
651 dealii::make_vectorized_array<Number>(0.5) * (uM * wM_n + uP * wP_n);
653 vector jump_value = uM - uP;
655 scalar lambda = calculate_lambda(wM_n, wP_n);
657 return (average_normal_flux + 0.5 * lambda * jump_value);
663 inline DEAL_II_ALWAYS_INLINE
665 calculate_lax_friedrichs_flux_linearized(vector
const & uM,
667 vector
const & delta_uM,
668 vector
const & delta_uP,
669 vector
const & normalM)
const
671 scalar uM_n = uM * normalM;
672 scalar uP_n = uP * normalM;
674 scalar delta_uM_n = delta_uM * normalM;
675 scalar delta_uP_n = delta_uP * normalM;
677 vector average_normal_flux =
678 dealii::make_vectorized_array<Number>(0.5) *
679 (uM * delta_uM_n + delta_uM * uM_n + uP * delta_uP_n + delta_uP * uP_n);
681 vector jump_value = delta_uM - delta_uP;
683 scalar lambda = calculate_lambda(uM_n, uP_n);
685 return (average_normal_flux + 0.5 * lambda * jump_value);
691 inline DEAL_II_ALWAYS_INLINE
693 calculate_upwind_flux(vector
const & uM,
695 scalar
const & average_normal_velocity)
const
697 vector average_velocity = 0.5 * (uM + uP);
699 vector jump_value = uM - uP;
701 return (average_normal_velocity * average_velocity +
702 data.upwind_factor * 0.5 * std::abs(average_normal_velocity) * jump_value);
708 inline DEAL_II_ALWAYS_INLINE
710 apply_outflow_bc(vector & flux, scalar
const & uM_n)
const
716 scalar outflow_indicator = dealii::make_vectorized_array<Number>(1.0);
718 for(
unsigned int v = 0; v < dealii::VectorizedArray<Number>::size(); ++v)
721 outflow_indicator[v] = 0.0;
725 flux = outflow_indicator * flux;
731 inline DEAL_II_ALWAYS_INLINE
733 calculate_upwind_flux_linearized(vector
const & uM,
735 vector
const & u_grid,
736 vector
const & delta_uM,
737 vector
const & delta_uP,
738 vector
const & normalM)
const
740 vector average_velocity = 0.5 * (uM + uP);
741 vector delta_average_velocity = 0.5 * (delta_uM + delta_uP);
743 scalar average_normal_velocity = average_velocity * normalM;
745 average_normal_velocity -= u_grid * normalM;
747 scalar delta_average_normal_velocity = delta_average_velocity * normalM;
749 vector jump_value = delta_uM - delta_uP;
751 return (average_normal_velocity * delta_average_velocity +
752 delta_average_normal_velocity * average_velocity +
753 data.upwind_factor * 0.5 * std::abs(average_normal_velocity) * jump_value);
764 inline DEAL_II_ALWAYS_INLINE
765 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>>
766 calculate_exterior_value_linearized(
767 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> & delta_uM,
768 unsigned int const q,
769 FaceIntegrator<dim, dim, Number> & integrator,
770 BoundaryTypeU
const & boundary_type)
const
773 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> delta_uP;
775 if(boundary_type == BoundaryTypeU::Dirichlet or boundary_type == BoundaryTypeU::DirichletCached)
777 if(data.type_dirichlet_bc == TypeDirichletBCs::Mirror)
779 delta_uP = -delta_uM;
781 else if(data.type_dirichlet_bc == TypeDirichletBCs::Direct)
791 "Type of imposition of Dirichlet BC's for convective term is not implemented."));
794 else if(boundary_type == BoundaryTypeU::Neumann)
798 else if(boundary_type == BoundaryTypeU::Symmetry)
800 dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> normalM =
801 integrator.get_normal_vector(q);
802 delta_uP = delta_uM - 2. * (delta_uM * normalM) * normalM;
807 dealii::ExcMessage(
"Boundary type of face is invalid or not implemented."));
819 std::shared_ptr<IntegratorCell> integrator_velocity;
820 std::shared_ptr<IntegratorFace> integrator_velocity_m;
821 std::shared_ptr<IntegratorFace> integrator_velocity_p;
823 std::shared_ptr<IntegratorCell> integrator_grid_velocity;
824 std::shared_ptr<IntegratorFace> integrator_grid_velocity_face;
860 typedef dealii::VectorizedArray<Number> scalar;
861 typedef dealii::Tensor<1, dim, dealii::VectorizedArray<Number>> vector;
862 typedef dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> tensor;
868 typedef typename Base::VectorType VectorType;
869 typedef typename Base::Range Range;
870 typedef typename Base::IntegratorCell IntegratorCell;
871 typedef typename Base::IntegratorFace IntegratorFace;
878 set_velocity_copy(VectorType
const & src)
const;
881 set_velocity_ptr(VectorType
const & src)
const;
883 dealii::LinearAlgebra::distributed::Vector<Number>
const &
884 get_velocity()
const;
887 initialize(dealii::MatrixFree<dim, Number>
const & matrix_free,
888 dealii::AffineConstraints<Number>
const & affine_constraints,
896 evaluate_nonlinear_operator(VectorType & dst, VectorType
const & src, Number
const time)
const;
899 evaluate_nonlinear_operator_add(VectorType & dst,
900 VectorType
const & src,
901 Number
const time)
const;
906 rhs(VectorType & dst)
const final;
909 rhs_add(VectorType & dst)
const final;
912 evaluate(VectorType & dst, VectorType
const & src)
const final;
915 evaluate_add(VectorType & dst, VectorType
const & src)
const final;
922 cell_loop_nonlinear_operator(dealii::MatrixFree<dim, Number>
const & matrix_free,
924 VectorType
const & src,
925 Range
const & cell_range)
const;
928 face_loop_nonlinear_operator(dealii::MatrixFree<dim, Number>
const & matrix_free,
930 VectorType
const & src,
931 Range
const & face_range)
const;
934 boundary_face_loop_nonlinear_operator(dealii::MatrixFree<dim, Number>
const & matrix_free,
936 VectorType
const & src,
937 Range
const & face_range)
const;
940 do_cell_integral_nonlinear_operator(IntegratorCell & integrator,
941 IntegratorCell & integrator_u_grid)
const;
944 do_face_integral_nonlinear_operator(IntegratorFace & integrator_m,
945 IntegratorFace & integrator_p,
946 IntegratorFace & integrator_grid_velocity)
const;
949 do_boundary_integral_nonlinear_operator(IntegratorFace & integrator,
950 IntegratorFace & integrator_grid_velocity,
951 dealii::types::boundary_id
const & boundary_id)
const;
959 reinit_cell_derived(IntegratorCell & integrator,
unsigned int const cell)
const final;
963 reinit_face_derived(IntegratorFace & integrator_m,
964 IntegratorFace & integrator_p,
965 unsigned int const face)
const final;
969 reinit_boundary_face_derived(IntegratorFace & integrator_m,
unsigned int const face)
const final;
973 reinit_face_cell_based_derived(IntegratorFace & integrator_m,
974 IntegratorFace & integrator_p,
975 unsigned int const cell,
976 unsigned int const face,
977 dealii::types::boundary_id
const boundary_id)
const final;
981 do_cell_integral(IntegratorCell & integrator)
const final;
985 do_face_integral(IntegratorFace & integrator_m, IntegratorFace & integrator_p)
const final;
989 do_face_int_integral(IntegratorFace & integrator_m, IntegratorFace & integrator_p)
const final;
998 do_face_int_integral_cell_based(IntegratorFace & integrator_m,
999 IntegratorFace & integrator_p)
const final;
1003 do_face_ext_integral(IntegratorFace & integrator_m, IntegratorFace & integrator_p)
const final;
1007 do_boundary_integral(IntegratorFace & integrator,
1008 OperatorType
const & operator_type,
1009 dealii::types::boundary_id
const & boundary_id)
const final;
1013 std::shared_ptr<Operators::ConvectiveKernel<dim, Number>> kernel;