ryujin 2.1.1 revision 9dcb748690310d6a540ebb8b066d1a0834fc7604
hyperbolic_system.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: MIT
3// Copyright (C) 2020 - 2023 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
10#include <discretization.h>
12#include <patterns_conversion.h>
13#include <simd.h>
14
15
16#include <deal.II/base/parameter_acceptor.h>
17#include <deal.II/base/tensor.h>
18
19#include <array>
20#include <functional>
21
22namespace ryujin
23{
24 namespace ShallowWater
25 {
35 class HyperbolicSystem final : public dealii::ParameterAcceptor
36 {
37 public:
41 static const std::string problem_name;
42
46 template <int dim>
47 static constexpr unsigned int problem_dimension = 1 + dim;
48
53 template <int dim, typename Number>
54 using state_type = dealii::Tensor<1, problem_dimension<dim>, Number>;
55
59 template <int dim>
60 static const std::array<std::string, problem_dimension<dim>>
62
66 template <int dim, typename Number>
68 dealii::Tensor<1, problem_dimension<dim>, Number>;
69
73 template <int dim>
74 static const std::array<std::string, problem_dimension<dim>>
76
80 template <int dim, typename Number>
81 using flux_type = dealii::
82 Tensor<1, problem_dimension<dim>, dealii::Tensor<1, dim, Number>>;
83
87 template <int dim, typename Number>
89 std::tuple<state_type<dim, Number>, Number>;
90
94 HyperbolicSystem(const std::string &subsection = "/HyperbolicSystem");
95
100
104 template <int dim>
105 static constexpr unsigned int n_precomputed_initial_values = 1;
106
110 template <int dim, typename Number>
112 std::array<Number, n_precomputed_initial_values<dim>>;
113
117 template <int dim>
118 static const std::array<std::string, n_precomputed_initial_values<dim>>
120
124 template <int dim>
125 static constexpr unsigned int n_precomputed_values = 1;
126
130 template <int dim, typename Number>
131 using precomputed_type = std::array<Number, n_precomputed_values<dim>>;
132
136 template <int dim>
137 static const std::array<std::string, n_precomputed_values<dim>>
139
143 static constexpr unsigned int n_precomputation_cycles = 1;
144
148 template <unsigned int cycle,
149 typename Number,
150 typename ScalarNumber = typename get_value_type<Number>::type,
151 int problem_dim,
152 typename MCV,
153 typename SPARSITY>
154 void
155 precomputation(MCV &precomputed_values,
157 const SPARSITY &sparsity_simd,
158 unsigned int i) const;
159
161
165
170 template <int problem_dim, typename Number>
171 static Number
172 water_depth(const dealii::Tensor<1, problem_dim, Number> &U);
173
180 template <int problem_dim, typename Number>
181 Number
182 water_depth_sharp(const dealii::Tensor<1, problem_dim, Number> &U) const;
183
190 template <int problem_dim, typename Number>
192 const dealii::Tensor<1, problem_dim, Number> &U) const;
193
194 template <int problem_dim, typename Number>
196 const dealii::Tensor<1, problem_dim, Number> &U) const;
197
204 template <typename Number>
205 Number filter_dry_water_depth(const Number &h) const;
206
211 template <int problem_dim, typename Number>
212 static dealii::Tensor<1, problem_dim - 1, Number>
213 momentum(const dealii::Tensor<1, problem_dim, Number> &U);
214
222 template <int problem_dim, typename Number>
223 Number
224 kinetic_energy(const dealii::Tensor<1, problem_dim, Number> &U) const;
225
233 template <int problem_dim, typename Number>
234 Number pressure(const dealii::Tensor<1, problem_dim, Number> &U) const;
235
243 template <int problem_dim, typename Number>
244 Number
245 speed_of_sound(const dealii::Tensor<1, problem_dim, Number> &U) const;
246
251 template <int problem_dim, typename Number>
253 const dealii::Tensor<1, problem_dim, Number> &U) const;
254
259 template <int problem_dim, typename Number>
260 dealii::Tensor<1, problem_dim, Number> mathematical_entropy_derivative(
261 const dealii::Tensor<1, problem_dim, Number> &U) const;
262
268 template <int problem_dim, typename Number>
269 bool is_admissible(const dealii::Tensor<1, problem_dim, Number> &U) const;
270
272
276
282 template <int component, int problem_dim, typename Number>
283 dealii::Tensor<1, problem_dim, Number> prescribe_riemann_characteristic(
284 const dealii::Tensor<1, problem_dim, Number> &U,
285 const dealii::Tensor<1, problem_dim, Number> &U_bar,
286 const dealii::Tensor<1, problem_dim - 1, Number> &normal) const;
287
306 template <int problem_dim, typename Number, typename Lambda>
307 dealii::Tensor<1, problem_dim, Number> apply_boundary_conditions(
308 dealii::types::boundary_id id,
309 dealii::Tensor<1, problem_dim, Number> U,
310 const dealii::Tensor<1, problem_dim - 1, Number> &normal,
311 Lambda get_dirichlet_data) const;
312
314
318
324 template <int problem_dim, typename Number>
325 dealii::Tensor<1, problem_dim, Number>
326 star_state(const dealii::Tensor<1, problem_dim, Number> &U,
327 const Number &Z_left,
328 const Number &Z_right) const;
329
339 template <typename ST,
340 int dim = ST::dimension - 1,
341 typename T = typename ST::value_type>
342 flux_type<dim, T> f(const ST &U) const;
343
361 template <typename MultiComponentVector,
362 typename MultiComponentVector2,
363 typename ST,
364 int dim = ST::dimension - 1,
365 typename T = typename ST::value_type>
367 flux_contribution(const MultiComponentVector &precomputed_values,
368 const MultiComponentVector2 &precomputed_initial_values,
369 const unsigned int i,
370 const ST &U_i) const;
371
372 template <typename MultiComponentVector,
373 typename MultiComponentVector2,
374 typename ST,
375 int dim = ST::dimension - 1,
376 typename T = typename ST::value_type>
378 flux_contribution(const MultiComponentVector &precomputed_values,
379 const MultiComponentVector2 &precomputed_initial_values,
380 const unsigned int *js,
381 const ST &U_j) const;
382
388 template <typename ST,
389 int dim = ST::dimension - 1,
390 typename T = typename ST::value_type>
391 flux_type<dim, T> flux(const std::tuple<ST, T> &prec_i,
392 const std::tuple<ST, T> &prec_j) const;
393
394 static constexpr bool have_high_order_flux = true;
395
400 template <typename ST,
401 int dim = ST::dimension - 1,
402 typename T = typename ST::value_type>
403 flux_type<dim, T> high_order_flux(const std::tuple<ST, T> &prec_i,
404 const std::tuple<ST, T> &prec_j) const;
405
406 static constexpr bool have_equilibrated_states = true;
407
413 template <typename ST,
414 int dim = ST::dimension - 1,
415 typename T = typename ST::value_type>
416 std::array<ST, 2>
417 equilibrated_states(const std::tuple<ST, T> &prec_i,
418 const std::tuple<ST, T> &prec_j) const;
419
421
425
426 static constexpr bool have_source_terms = true;
427
431 template <typename MultiComponentVector,
432 typename ST,
433 int dim = ST::dimension - 1,
434 typename T = typename ST::value_type>
435 ST low_order_nodal_source(const MultiComponentVector &precomputed_values,
436 const unsigned int i,
437 const ST &U_i) const;
438
442 template <typename MultiComponentVector,
443 typename ST,
444 int dim = ST::dimension - 1,
445 typename T = typename ST::value_type>
446 ST high_order_nodal_source(const MultiComponentVector &precomputed_values,
447 const unsigned int i,
448 const ST &U_i) const;
449
455 template <typename ST,
456 int dim = ST::dimension - 1,
457 typename T = typename ST::value_type>
458 ST low_order_stencil_source(const std::tuple<ST, T> &prec_i,
459 const std::tuple<ST, T> &prec_j,
460 const T &d_ij,
461 const dealii::Tensor<1, dim, T> &c_ij) const;
462
463
468 template <typename ST,
469 int dim = ST::dimension - 1,
470 typename T = typename ST::value_type>
471 ST high_order_stencil_source(const std::tuple<ST, T> &prec_i,
472 const std::tuple<ST, T> &prec_j,
473 const T &d_ij,
474 const dealii::Tensor<1, dim, T> &c_ij) const;
475
476
482 template <typename ST,
483 int dim = ST::dimension - 1,
484 typename T = typename ST::value_type>
485 ST
486 affine_shift_stencil_source(const std::tuple<ST, T> &prec_i,
487 const std::tuple<ST, T> &prec_j,
488 const T &d_ij,
489 const dealii::Tensor<1, dim, T> &c_ij) const;
490
491
493
498
499 /*
500 * Given a state vector associated with @ref dim2 spatial dimensions
501 * return an "expanded" version of the state vector associated with
502 * @ref dim1 spatial dimensions where the momentum vector is projected
503 * onto the first @ref dim2 unit directions of the @ref dim1
504 * dimensional euclidean space.
505 *
506 * @precondition dim has to be larger or equal than dim2.
507 */
508 template <int dim,
509 typename ST,
510 typename T = typename ST::value_type,
511 int dim2 = ST::dimension - 1,
512 typename = typename std::enable_if<(dim >= dim2)>::type>
513 state_type<dim, T> expand_state(const ST &state) const;
514
515 /*
516 * Given a primitive state [h, u_1, ..., u_d] return a conserved
517 * state
518 */
519 template <int problem_dim, typename Number>
520 dealii::Tensor<1, problem_dim, Number> from_primitive_state(
521 const dealii::Tensor<1, problem_dim, Number> &primitive_state) const;
522
523 /*
524 * Given a conserved state return a primitive state [h, u_1, ..., u_d]
525 */
526 template <int problem_dim, typename Number>
527 dealii::Tensor<1, problem_dim, Number> to_primitive_state(
528 const dealii::Tensor<1, problem_dim, Number> &state) const;
529
530 /*
531 * Transform the current state according to a given operator @ref
532 * momentum_transform acting on a @p dim dimensional momentum (or
533 * velocity) vector.
534 */
535 template <int problem_dim, typename Number, typename Lambda>
536 dealii::Tensor<1, problem_dim, Number> apply_galilei_transform(
537 const dealii::Tensor<1, problem_dim, Number> &state,
538 const Lambda &lambda) const;
539
541
545
546 ACCESSOR_READ_ONLY(gravity)
547
548 ACCESSOR_READ_ONLY(mannings)
549
550 ACCESSOR_READ_ONLY(reference_water_depth)
551
552 ACCESSOR_READ_ONLY(dry_state_relaxation_sharp)
553
554 ACCESSOR_READ_ONLY(dry_state_relaxation_mollified)
555
556 private:
557 double gravity_;
558
559 double mannings_;
560
561 double reference_water_depth_;
562
563 double dry_state_relaxation_sharp_;
564
565 double dry_state_relaxation_mollified_;
566
568 };
569
570
571 /* Inline definitions */
572
573
574 template <unsigned int cycle,
575 typename Number,
576 typename ScalarNumber,
577 int problem_dim,
578 typename MCV,
579 typename SPARSITY>
580 DEAL_II_ALWAYS_INLINE inline void HyperbolicSystem::precomputation(
581 MCV &precomputed_values,
583 const SPARSITY & /*sparsity_simd*/,
584 unsigned int i) const
585 {
586 static_assert(cycle == 0, "internal error");
587
588 constexpr int dim = problem_dim - 1;
589
590 const auto U_i = U.template get_tensor<Number>(i);
591
593 precomputed_values.template write_tensor<Number>(prec_i, i);
594 }
595
596
597 template <int problem_dim, typename Number>
598 DEAL_II_ALWAYS_INLINE inline Number HyperbolicSystem::water_depth(
599 const dealii::Tensor<1, problem_dim, Number> &U)
600 {
601 return U[0];
602 }
603
604
605 template <int problem_dim, typename Number>
606 DEAL_II_ALWAYS_INLINE inline Number
608 const dealii::Tensor<1, problem_dim, Number> &U) const
609 {
610 using ScalarNumber = typename get_value_type<Number>::type;
611 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
612 const Number h_cutoff_mollified =
613 Number(reference_water_depth_ * dry_state_relaxation_mollified_) *
614 eps;
615
616 const Number h = positive_part(water_depth(U));
617 const Number h_max = std::max(h, h_cutoff_mollified);
618 const Number denom = h * h + h_max * h_max;
619 return ScalarNumber(2.) * h / denom;
620 }
621
622
623 template <int problem_dim, typename Number>
624 DEAL_II_ALWAYS_INLINE inline Number HyperbolicSystem::water_depth_sharp(
625 const dealii::Tensor<1, problem_dim, Number> &U) const
626 {
627 using ScalarNumber = typename get_value_type<Number>::type;
628 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
629 const Number h_cutoff_sharp =
630 Number(reference_water_depth_ * dry_state_relaxation_sharp_) * eps;
631
632 const Number h = water_depth(U);
633 const Number h_max = std::max(h, h_cutoff_sharp);
634 return h_max;
635 }
636
637
638 template <int problem_dim, typename Number>
639 DEAL_II_ALWAYS_INLINE inline Number
641 const dealii::Tensor<1, problem_dim, Number> &U) const
642 {
643 using ScalarNumber = typename get_value_type<Number>::type;
644 return ScalarNumber(1.) / water_depth_sharp(U);
645 }
646
647
648 template <typename Number>
649 DEAL_II_ALWAYS_INLINE inline Number
651 {
652 using ScalarNumber = typename get_value_type<Number>::type;
653 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
654 const Number h_cutoff_mollified =
655 Number(reference_water_depth_ * dry_state_relaxation_mollified_) *
656 eps;
657
658 return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
659 std::abs(h), h_cutoff_mollified, Number(0.), h);
660 }
661
662
663 template <int problem_dim, typename Number>
664 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, problem_dim - 1, Number>
665 HyperbolicSystem::momentum(const dealii::Tensor<1, problem_dim, Number> &U)
666 {
667 constexpr int dim = problem_dim - 1;
668
669 dealii::Tensor<1, dim, Number> result;
670
671 for (unsigned int i = 0; i < dim; ++i)
672 result[i] = U[1 + i];
673 return result;
674 }
675
676
677 template <int problem_dim, typename Number>
678 DEAL_II_ALWAYS_INLINE inline Number HyperbolicSystem::kinetic_energy(
679 const dealii::Tensor<1, problem_dim, Number> &U) const
680 {
681 /* KE = 1/2 h |v|^2 */
682
683 using ScalarNumber = typename get_value_type<Number>::type;
684
685 const auto h = water_depth(U);
686 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
687
688 return ScalarNumber(0.5) * h * vel.norm_square();
689 }
690
691
692 template <int problem_dim, typename Number>
693 DEAL_II_ALWAYS_INLINE inline Number HyperbolicSystem::pressure(
694 const dealii::Tensor<1, problem_dim, Number> &U) const
695 {
696 /* p = 1/2 g h^2 */
697
698 using ScalarNumber = typename get_value_type<Number>::type;
699
700 const Number h_sqd = U[0] * U[0];
701
702 return ScalarNumber(0.5 * gravity_) * h_sqd;
703 }
704
705
706 template <int problem_dim, typename Number>
707 DEAL_II_ALWAYS_INLINE inline Number HyperbolicSystem::speed_of_sound(
708 const dealii::Tensor<1, problem_dim, Number> &U) const
709 {
710 using ScalarNumber = typename get_value_type<Number>::type;
711
712 /* c^2 = g * h */
713 return std::sqrt(ScalarNumber(gravity_) * U[0]);
714 }
715
716
717 template <int problem_dim, typename Number>
718 DEAL_II_ALWAYS_INLINE inline Number HyperbolicSystem::mathematical_entropy(
719 const dealii::Tensor<1, problem_dim, Number> &U) const
720 {
721 const auto p = pressure(U);
722 const auto k_e = kinetic_energy(U);
723 return p + k_e;
724 }
725
726
727 template <int problem_dim, typename Number>
728 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, problem_dim, Number>
730 const dealii::Tensor<1, problem_dim, Number> &U) const
731 {
732 /*
733 * With
734 * eta = 1/2 g h^2 + 1/2 |m|^2 / h
735 *
736 * we get
737 *
738 * eta' = (g h - 1/2 |vel|^2, vel)
739 *
740 * where vel = m / h
741 */
742
743 constexpr int dim = problem_dim - 1;
744 using ScalarNumber = typename get_value_type<Number>::type;
745
746 dealii::Tensor<1, problem_dim, Number> result;
747
748 const Number &h = U[0];
749 const auto vel = momentum(U) * inverse_water_depth_sharp(U);
750
751 // water depth component
752 result[0] =
753 ScalarNumber(gravity_) * h - ScalarNumber(0.5) * vel.norm_square();
754
755 // momentum components
756 for (unsigned int i = 0; i < dim; ++i) {
757 result[1 + i] = vel[i];
758 }
759
760 return result;
761 }
762
763
764 template <int problem_dim, typename T>
765 DEAL_II_ALWAYS_INLINE inline bool HyperbolicSystem::is_admissible(
766 const dealii::Tensor<1, problem_dim, T> &U) const
767 {
768 const auto h = filter_dry_water_depth(water_depth(U));
769
770 constexpr auto gte = dealii::SIMDComparison::greater_than_or_equal;
771 const auto test =
772 dealii::compare_and_apply_mask<gte>(h, T(0.), T(0.), T(-1.));
773
774#ifdef DEBUG_OUTPUT
775 if (!(test == T(0.))) {
776 std::cout << std::fixed << std::setprecision(16);
777 std::cout << "Bounds violation: Negative state [h] detected!\n";
778 std::cout << "\t\th: " << h << "\n" << std::endl;
779 }
780#endif
781
782 return (test == T(0.));
783 }
784
785
786 template <int component, int problem_dim, typename Number>
787 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, problem_dim, Number>
789 const dealii::Tensor<1, problem_dim, Number> &U,
790 const dealii::Tensor<1, problem_dim, Number> &U_bar,
791 const dealii::Tensor<1, problem_dim - 1, Number> &normal) const
792 {
793
794 constexpr int dim = problem_dim - 1;
795
796 /* Note that U_bar are the dirichlet values that are prescribed */
797 static_assert(component == 1 || component == 2,
798 "component has to be 1 or 2");
799
800 using ScalarNumber = typename get_value_type<Number>::type;
801
802 const auto m = momentum(U);
803 const auto a = speed_of_sound(U);
804 const auto vn = m * normal * inverse_water_depth_sharp(U);
805
806 const auto m_bar = momentum(U_bar);
807 const auto a_bar = speed_of_sound(U_bar);
808 const auto vn_bar = m_bar * normal * inverse_water_depth_sharp(U_bar);
809
810 /* First Riemann characteristic: v * n - 2 * a */
811
812 const auto R_1 = component == 1 ? vn_bar - ScalarNumber(2.) * a_bar
813 : vn - ScalarNumber(2.) * a;
814
815 /* Second Riemann characteristic: v * n + 2 * a */
816
817 const auto R_2 = component == 2 ? vn_bar + ScalarNumber(2.) * a_bar
818 : vn + ScalarNumber(2.) * a;
819
820 const auto vperp = m * inverse_water_depth_sharp(U) - vn * normal;
821
822 const auto vn_new = ScalarNumber(0.5) * (R_1 + R_2);
823
824 const auto h_new =
825 ryujin::fixed_power<2>((R_2 - R_1) / ScalarNumber(4.)) /
826 ScalarNumber(gravity_);
827
829
830 U_new[0] = h_new;
831
832 for (unsigned int d = 0; d < dim; ++d) {
833 U_new[1 + d] = h_new * (vn_new * normal + vperp)[d];
834 }
835
836 return U_new;
837 }
838
839
840 template <int problem_dim, typename Number, typename Lambda>
841 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, problem_dim, Number>
843 dealii::types::boundary_id id,
844 dealii::Tensor<1, problem_dim, Number> U,
845 const dealii::Tensor<1, problem_dim - 1, Number> &normal,
846 Lambda get_dirichlet_data) const
847 {
848 constexpr auto dim = problem_dim - 1;
849
850 if (id == Boundary::dirichlet) {
851 U = get_dirichlet_data();
852
853 } else if (id == Boundary::slip) {
854 auto m = momentum(U);
855 m -= 1. * (m * normal) * normal;
856 for (unsigned int k = 0; k < dim; ++k)
857 U[k + 1] = m[k];
858
859 } else if (id == Boundary::no_slip) {
860 for (unsigned int k = 0; k < dim; ++k)
861 U[k + 1] = Number(0.);
862
863 } else if (id == Boundary::dynamic) {
864 /*
865 * On dynamic boundary conditions, we distinguish four cases:
866 *
867 * - supersonic inflow: prescribe full state
868 * - subsonic inflow:
869 * decompose into Riemann invariants and leave R_2
870 * characteristic untouched.
871 * - supersonic outflow: do nothing
872 * - subsonic outflow:
873 * decompose into Riemann invariants and prescribe incoming
874 * R_1 characteristic.
875 */
876 const auto m = momentum(U);
877 const auto h_inverse = inverse_water_depth_sharp(U);
878 const auto a = speed_of_sound(U);
879 const auto vn = m * normal * h_inverse;
880
881 /* Supersonic inflow: */
882 if (vn < -a) {
883 U = get_dirichlet_data();
884 }
885
886 /* Subsonic inflow: */
887 if (vn >= -a && vn <= 0.) {
888 const auto U_dirichlet = get_dirichlet_data();
889 U = prescribe_riemann_characteristic<2>(U_dirichlet, U, normal);
890 }
891
892 /* Subsonic outflow: */
893 if (vn > 0. && vn <= a) {
894 const auto U_dirichlet = get_dirichlet_data();
895 U = prescribe_riemann_characteristic<1>(U, U_dirichlet, normal);
896 }
897
898 /* Supersonic outflow: do nothing, i.e., keep U as is */
899 }
900
901 return U;
902 }
903
904
905 template <int problem_dim, typename Number>
906 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, problem_dim, Number>
908 const dealii::Tensor<1, problem_dim, Number> &U_left,
909 const Number &Z_left,
910 const Number &Z_right) const
911 {
912 const Number Z_max = std::max(Z_left, Z_right);
913 const Number h = water_depth(U_left);
914 const Number H_star = std::max(Number(0.), h + Z_left - Z_max);
915
916 return U_left * H_star * inverse_water_depth_mollified(U_left);
917 }
918
919
920 template <typename ST, int dim, typename T>
921 DEAL_II_ALWAYS_INLINE inline auto HyperbolicSystem::f(const ST &U) const
923 {
924 const T h_inverse = inverse_water_depth_sharp(U);
925 const auto m = momentum(U);
926 const auto p = pressure(U);
927
928 flux_type<dim, T> result;
929
930 result[0] = m;
931 for (unsigned int i = 0; i < dim; ++i) {
932 result[1 + i] = m * (m[i] * h_inverse);
933 result[1 + i][i] += p;
934 }
935 return result;
936 }
937
938
939 template <typename MCV, typename MICV, typename ST, int dim, typename T>
940 DEAL_II_ALWAYS_INLINE inline auto
941 HyperbolicSystem::flux_contribution(const MCV & /*precomputed_values*/,
942 const MICV &precomputed_initial_values,
943 const unsigned int i,
944 const ST &U_i) const
946 {
947 const auto Z_i = precomputed_initial_values.template get_tensor<T>(i)[0];
948 return {U_i, Z_i};
949 }
950
951
952 template <typename MCV, typename MICV, typename ST, int dim, typename T>
953 DEAL_II_ALWAYS_INLINE inline auto
954 HyperbolicSystem::flux_contribution(const MCV & /*precomputed_values*/,
955 const MICV &precomputed_initial_values,
956 const unsigned int *js,
957 const ST &U_j) const
959 {
960 const auto Z_j = precomputed_initial_values.template get_tensor<T>(js)[0];
961 return {U_j, Z_j};
962 }
963
964
965 template <typename ST, int dim, typename T>
966 DEAL_II_ALWAYS_INLINE inline auto
967 HyperbolicSystem::flux(const std::tuple<ST, T> &prec_i,
968 const std::tuple<ST, T> &prec_j) const
970 {
971 const auto &[U_star_ij, U_star_ji] = equilibrated_states(prec_i, prec_j);
972
973 const auto f_i = f(U_star_ij);
974 const auto f_j = f(U_star_ji);
975
976 return -add(f_i, f_j);
977 }
978
979
980 template <typename ST, int dim, typename T>
981 DEAL_II_ALWAYS_INLINE inline auto
982 HyperbolicSystem::high_order_flux(const std::tuple<ST, T> &prec_i,
983 const std::tuple<ST, T> &prec_j) const
985 {
986 const auto &[U_i, Z_i] = prec_i;
987 const auto &[U_j, Z_j] = prec_j;
988
989 const auto f_i = f(U_i);
990 const auto f_j = f(U_j);
991
992 return -add(f_i, f_j);
993 }
994
995
996 template <typename ST, int dim, typename T>
997 DEAL_II_ALWAYS_INLINE inline std::array<ST, 2>
998 HyperbolicSystem::equilibrated_states(const std::tuple<ST, T> &prec_i,
999 const std::tuple<ST, T> &prec_j) const
1000 {
1001 const auto &[U_i, Z_i] = prec_i;
1002 const auto &[U_j, Z_j] = prec_j;
1003 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1004 const auto U_star_ji = star_state(U_j, Z_j, Z_i);
1005 return {U_star_ij, U_star_ji};
1006 }
1007
1008
1009 template <typename MultiComponentVector, typename ST, int dim, typename T>
1011 const MultiComponentVector & /*precomputed_values*/,
1012 const unsigned int /*i*/,
1013 const ST & /*U_i*/) const
1014 {
1015 // FIXME
1016 return ST();
1017 }
1018
1019
1020 template <typename MultiComponentVector, typename ST, int dim, typename T>
1022 const MultiComponentVector & /*precomputed_values*/,
1023 const unsigned int /*i*/,
1024 const ST & /*U_i*/) const
1025 {
1026 // FIXME
1027 return ST();
1028 }
1029
1030
1031 template <typename ST, int dim, typename T>
1033 const std::tuple<ST, T> &prec_i,
1034 const std::tuple<ST, T> &prec_j,
1035 const T &,
1036 const dealii::Tensor<1, dim, T> &c_ij) const
1037 {
1038 using ScalarNumber = typename get_value_type<T>::type;
1039
1040 const auto &[U_i, Z_i] = prec_i;
1041 const auto H_i = water_depth(U_i);
1042 const auto &[U_j, Z_j] = prec_j;
1043 const auto U_star_ij = star_state(U_i, Z_i, Z_j);
1044 const auto H_star_ij = water_depth(U_star_ij);
1045
1046 const auto factor =
1047 ScalarNumber(gravity_) * (H_star_ij * H_star_ij - H_i * H_i);
1048 ST result;
1049 for (unsigned int d = 1; d < dim + 1; ++d)
1050 result[d] = factor * c_ij[d - 1];
1051 return result;
1052 }
1053
1054
1055 template <typename ST, int dim, typename T>
1057 const std::tuple<ST, T> &prec_i,
1058 const std::tuple<ST, T> &prec_j,
1059 const T &,
1060 const dealii::Tensor<1, dim, T> &c_ij) const
1061 {
1062 using ScalarNumber = typename get_value_type<T>::type;
1063
1064 const auto &[U_i, Z_i] = prec_i;
1065 const auto &[U_j, Z_j] = prec_j;
1066 const auto H_i = water_depth(U_i);
1067 const auto H_j = water_depth(U_j);
1068
1069 const auto factor = -ScalarNumber(gravity_) * H_i * (Z_j - Z_i) +
1070 ScalarNumber(0.5) * ScalarNumber(gravity_) *
1071 (H_j - H_i) * (H_j - H_i);
1072
1073 ST result;
1074 for (unsigned int d = 1; d < dim + 1; ++d)
1075 result[d] = factor * c_ij[d - 1];
1076 return result;
1077 }
1078
1079
1080 template <typename ST, int dim, typename T>
1082 const std::tuple<ST, T> &prec_i,
1083 const std::tuple<ST, T> &prec_j,
1084 const T &d_ij,
1085 const dealii::Tensor<1, dim, T> &c_ij) const
1086 {
1087 using ScalarNumber = typename get_value_type<T>::type;
1088
1089 const auto &[U_star_ij, U_star_ji] = equilibrated_states(prec_i, prec_j);
1090 const auto f_star_ij = f(U_star_ij);
1091
1092 return -ScalarNumber(2.) * d_ij * U_star_ij -
1093 ScalarNumber(2.) * contract(f_star_ij, c_ij);
1094 }
1095
1096
1097 template <int dim, typename ST, typename T, int dim2, typename>
1098 auto HyperbolicSystem::expand_state(const ST &state) const
1100 {
1101 state_type<dim, T> result;
1102 result[0] = state[0];
1103 for (unsigned int i = 1; i < dim2 + 1; ++i)
1104 result[i] = state[i];
1105
1106 return result;
1107 }
1108
1109
1110 template <int problem_dim, typename Number>
1111 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, problem_dim, Number>
1113 const dealii::Tensor<1, problem_dim, Number> &primitive_state) const
1114 {
1115 const auto &h = primitive_state[0];
1116
1117 auto state = primitive_state;
1118 /* Fix up momentum: */
1119 for (unsigned int i = 1; i < problem_dim; ++i)
1120 state[i] *= h;
1121
1122 return state;
1123 }
1124
1125
1126 template <int problem_dim, typename Number>
1127 DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, problem_dim, Number>
1129 const dealii::Tensor<1, problem_dim, Number> &state) const
1130 {
1131 const auto h_inverse = inverse_water_depth_sharp(state);
1132
1133 auto primitive_state = state;
1134 /* Fix up velocity: */
1135 for (unsigned int i = 1; i < problem_dim; ++i)
1136 primitive_state[i] *= h_inverse;
1137
1138 return primitive_state;
1139 }
1140
1141
1142 template <int problem_dim, typename Number, typename Lambda>
1143 dealii::Tensor<1, problem_dim, Number>
1145 const dealii::Tensor<1, problem_dim, Number> &state,
1146 const Lambda &lambda) const
1147 {
1148 constexpr auto dim = problem_dim - 1;
1149
1150 auto result = state;
1151 auto M = lambda(momentum(state));
1152 for (unsigned int d = 0; d < dim; ++d)
1153 result[1 + d] = M[d];
1154 return result;
1155 }
1156
1157 } // namespace ShallowWater
1158} // namespace ryujin
static Number water_depth(const dealii::Tensor< 1, problem_dim, Number > &U)
ST low_order_nodal_source(const MultiComponentVector &precomputed_values, const unsigned int i, const ST &U_i) const
ST high_order_stencil_source(const std::tuple< ST, T > &prec_i, const std::tuple< ST, T > &prec_j, const T &d_ij, const dealii::Tensor< 1, dim, T > &c_ij) const
ST high_order_nodal_source(const MultiComponentVector &precomputed_values, const unsigned int i, const ST &U_i) const
dealii::Tensor< 1, problem_dimension< dim >, Number > primitive_state_type
static constexpr unsigned int n_precomputed_initial_values
dealii::Tensor< 1, problem_dim, Number > from_primitive_state(const dealii::Tensor< 1, problem_dim, Number > &primitive_state) const
static dealii::Tensor< 1, problem_dim - 1, Number > momentum(const dealii::Tensor< 1, problem_dim, Number > &U)
dealii::Tensor< 1, problem_dim, Number > star_state(const dealii::Tensor< 1, problem_dim, Number > &U, const Number &Z_left, const Number &Z_right) const
static constexpr unsigned int problem_dimension
void precomputation(MCV &precomputed_values, const MultiComponentVector< ScalarNumber, problem_dim > &U, const SPARSITY &sparsity_simd, unsigned int i) const
static constexpr unsigned int n_precomputed_values
static constexpr bool have_equilibrated_states
static const std::array< std::string, n_precomputed_initial_values< dim > > precomputed_initial_names
Number inverse_water_depth_mollified(const dealii::Tensor< 1, problem_dim, Number > &U) const
HyperbolicSystem(const std::string &subsection="/HyperbolicSystem")
static const std::array< std::string, problem_dimension< dim > > primitive_component_names
ST low_order_stencil_source(const std::tuple< ST, T > &prec_i, const std::tuple< ST, T > &prec_j, const T &d_ij, const dealii::Tensor< 1, dim, T > &c_ij) const
Number mathematical_entropy(const dealii::Tensor< 1, problem_dim, Number > &U) const
std::array< ST, 2 > equilibrated_states(const std::tuple< ST, T > &prec_i, const std::tuple< ST, T > &prec_j) const
dealii::Tensor< 1, problem_dim, Number > apply_boundary_conditions(dealii::types::boundary_id id, dealii::Tensor< 1, problem_dim, Number > U, const dealii::Tensor< 1, problem_dim - 1, Number > &normal, Lambda get_dirichlet_data) const
Number water_depth_sharp(const dealii::Tensor< 1, problem_dim, Number > &U) const
flux_contribution_type< dim, T > flux_contribution(const MultiComponentVector &precomputed_values, const MultiComponentVector2 &precomputed_initial_values, const unsigned int *js, const ST &U_j) const
dealii::Tensor< 1, problem_dim, Number > prescribe_riemann_characteristic(const dealii::Tensor< 1, problem_dim, Number > &U, const dealii::Tensor< 1, problem_dim, Number > &U_bar, const dealii::Tensor< 1, problem_dim - 1, Number > &normal) const
flux_type< dim, T > f(const ST &U) const
flux_type< dim, T > flux(const std::tuple< ST, T > &prec_i, const std::tuple< ST, T > &prec_j) const
static constexpr unsigned int n_precomputation_cycles
flux_type< dim, T > high_order_flux(const std::tuple< ST, T > &prec_i, const std::tuple< ST, T > &prec_j) const
Number speed_of_sound(const dealii::Tensor< 1, problem_dim, Number > &U) const
dealii::Tensor< 1, problem_dimension< dim >, dealii::Tensor< 1, dim, Number > > flux_type
Number pressure(const dealii::Tensor< 1, problem_dim, Number > &U) const
std::tuple< state_type< dim, Number >, Number > flux_contribution_type
flux_contribution_type< dim, T > flux_contribution(const MultiComponentVector &precomputed_values, const MultiComponentVector2 &precomputed_initial_values, const unsigned int i, const ST &U_i) const
Number kinetic_energy(const dealii::Tensor< 1, problem_dim, Number > &U) const
dealii::Tensor< 1, problem_dim, Number > to_primitive_state(const dealii::Tensor< 1, problem_dim, Number > &state) const
static const std::array< std::string, n_precomputed_values< dim > > precomputed_names
state_type< dim, T > expand_state(const ST &state) const
ST affine_shift_stencil_source(const std::tuple< ST, T > &prec_i, const std::tuple< ST, T > &prec_j, const T &d_ij, const dealii::Tensor< 1, dim, T > &c_ij) const
bool is_admissible(const dealii::Tensor< 1, problem_dim, Number > &U) const
static const std::array< std::string, problem_dimension< dim > > component_names
dealii::Tensor< 1, problem_dim, Number > mathematical_entropy_derivative(const dealii::Tensor< 1, problem_dim, Number > &U) const
std::array< Number, n_precomputed_values< dim > > precomputed_type
std::array< Number, n_precomputed_initial_values< dim > > precomputed_initial_type
dealii::Tensor< 1, problem_dim, Number > apply_galilei_transform(const dealii::Tensor< 1, problem_dim, Number > &state, const Lambda &lambda) const
dealii::Tensor< 1, problem_dimension< dim >, Number > state_type
Number inverse_water_depth_sharp(const dealii::Tensor< 1, problem_dim, Number > &U) const
Number filter_dry_water_depth(const Number &h) const
#define ACCESSOR_READ_ONLY(member)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
Definition: simd.h:111
DEAL_II_ALWAYS_INLINE FT add(const FT &flux_left_ij, const FT &flux_right_ij)
DEAL_II_ALWAYS_INLINE dealii::Tensor< 1, problem_dim, T > contract(const FT &flux_ij, const TT &c_ij)