[test.inviscid-compressible-flow] add another 1d euler test case

dune/gdt/test/instationary-eocstudies/base.hh

@@ -299,7 +299,7 @@ protected:
     return time_points;
   }
 
-  const double T_end_;
+  double T_end_;
   const size_t num_refinements_;
   const size_t num_additional_refinements_for_reference_;
   const std::function<void(const DiscreteBochnerFunction<V, GV, m>&, const std::string&)> visualize_;

dune/gdt/test/inviscid-compressible-flow/inviscid_compressible_flow__euler_1d__explicit__fv.cc

@@ -16,6 +16,8 @@
 #include <dune/xt/grid/grids.hh>
 #include <dune/xt/grid/walker.hh>
 #include <dune/xt/functions/base/sliced.hh>
+#include <dune/xt/functions/interfaces/element-functions.hh>
+#include <dune/xt/functions/lambda/function.hh>
 
 #include <dune/gdt/local/integrands/identity.hh>
 #include <dune/gdt/test/instationary-eocstudies/hyperbolic-fv.hh>
@@ -40,6 +42,7 @@ protected:
   using BaseType::m;
   using typename BaseType::E;
   using typename BaseType::F;
+  using typename BaseType::R;
   using typename BaseType::DF;
   using typename BaseType::GV;
   using typename BaseType::GP;
@@ -78,7 +81,7 @@ public:
 
   std::vector<std::string> quantities() const override final
   {
-    return {"rel mass conserv   error"};
+    return {"rel mass conserv   error"}; // <- This is on purpose, see column header formatting in ConvergenceStudy.
   }
 
   virtual std::map<std::string, std::map<std::string, double>>
@@ -124,8 +127,35 @@ protected:
 
   DF make_initial_values(const S& space) override final
   {
-    return Problem::access().template make_initial_values<V>(space);
-  }
+    if (boundary_treatment == "inflow_from_the_left_by_heuristic_euler_treatment_impermeable_wall_right") {
+      const auto& euler_tools = Problem::access().euler_tools;
+      return interpolate<V>(0,
+                            [&](const auto& /*xx*/, const auto& /*param*/) {
+                              FieldVector<R, m> primitive_variables(0.);
+                              // density
+                              primitive_variables[0] = 0.5;
+                              // velocity
+                              for (size_t ii = 0; ii < d; ++ii)
+                                primitive_variables[1 + ii] = 0.;
+                              // pressure
+                              primitive_variables[m - 1] = 0.4;
+                              return euler_tools.to_conservative(primitive_variables);
+                            },
+                            space);
+    } else
+      return Problem::access().template make_initial_values<V>(space);
+  } // ... make_initial_values(...)
+
+  double estimate_dt(const S& space) override final
+  {
+    if (boundary_treatment == "inflow_from_the_left_by_heuristic_euler_treatment_impermeable_wall_right")
+      return estimate_dt_for_hyperbolic_system(space.grid_view(),
+                                               make_initial_values(space),
+                                               flux(),
+                                               /*boundary_data_range=*/{{0.5, 0., 0.4}, {1.5, 0.5, 0.4}});
+    else
+      return estimate_dt_for_hyperbolic_system(space.grid_view(), make_initial_values(space), flux());
+  } // ... estimate_dt(...)
 
   GP make_initial_grid() override final
   {
@@ -146,17 +176,18 @@ protected:
       boundary_info = std::make_unique<XT::Grid::NormalBasedBoundaryInfo<XT::Grid::extract_intersection_t<GV>>>();
       boundary_info->register_new_normal({-1}, new XT::Grid::ImpermeableBoundary());
       boundary_info->register_new_normal({1}, new XT::Grid::ImpermeableBoundary());
-      XT::Grid::ApplyOn::CustomBoundaryIntersections<GV> impermeable_wall_filter(*boundary_info,
-                                                                                 new XT::Grid::ImpermeableBoundary());
+      const XT::Grid::ApplyOn::CustomBoundaryIntersections<GV> impermeable_wall_filter(
+          *boundary_info, new XT::Grid::ImpermeableBoundary());
       auto op = std::make_unique<AdvectionFvOperator<GV, V, m>>(space.grid_view(),
                                                                 *self.numerical_flux_,
                                                                 space,
                                                                 space,
                                                                 /*periodicity_restriction=*/impermeable_wall_filter);
       // the actual handling of impermeable walls
-      op->append([&](const auto& u,
+      op->append(/*numerical_boundary_flux=*/[&](
+                     const auto& u,
                      const auto& n,
-                     const auto& /*param*/ = {}) { return euler_tools.flux_at_impermeable_walls(u, n); },
+                     const auto& /*param*/) { return euler_tools.flux_at_impermeable_walls(u, n); },
                  {},
                  impermeable_wall_filter);
       return op;
@@ -164,8 +195,8 @@ protected:
       boundary_info = std::make_unique<XT::Grid::NormalBasedBoundaryInfo<XT::Grid::extract_intersection_t<GV>>>();
       boundary_info->register_new_normal({-1}, new XT::Grid::ImpermeableBoundary());
       boundary_info->register_new_normal({1}, new XT::Grid::ImpermeableBoundary());
-      XT::Grid::ApplyOn::CustomBoundaryIntersections<GV> impermeable_wall_filter(*boundary_info,
-                                                                                 new XT::Grid::ImpermeableBoundary());
+      const XT::Grid::ApplyOn::CustomBoundaryIntersections<GV> impermeable_wall_filter(
+          *boundary_info, new XT::Grid::ImpermeableBoundary());
       auto op = std::make_unique<AdvectionFvOperator<GV, V, m>>(space.grid_view(),
                                                                 *self.numerical_flux_,
                                                                 space,
@@ -173,11 +204,104 @@ protected:
                                                                 /*periodicity_restriction=*/impermeable_wall_filter);
       // the actual handling of impermeable walls, see [DF2015, p. 415, (8.66 - 8.67)]
       op->append(
+          /*boundary_extrapolation=*/
           [&](const auto& intersection,
               const auto& xx_in_reference_intersection_coordinates,
-              const auto& /*f*/,
+              const auto& /*flux*/,
               const auto& u,
-              const auto& /*mu*/ = {}) {
+              const auto& /*param*/) {
+            const auto normal = intersection.unitOuterNormal(xx_in_reference_intersection_coordinates);
+            const auto rho = euler_tools.density_from_conservative(u);
+            auto velocity = euler_tools.velocity_from_conservative(u);
+            velocity -= normal * 2. * (velocity * normal);
+            const auto pressure = euler_tools.pressure_from_conservative(u);
+            return euler_tools.to_conservative(XT::Common::hstack(rho, velocity, pressure));
+          },
+          {},
+          impermeable_wall_filter);
+      return op;
+    } else if (boundary_treatment == "inflow_from_the_left_by_heuristic_euler_treatment_impermeable_wall_right") {
+      periodic_density_variation_ =
+          std::unique_ptr<XT::Functions::LambdaFunction<d, m>>(new XT::Functions::LambdaFunction<d, m>(
+              /*order=*/0,
+              [&](const auto& /*xx*/, const auto& param) {
+                FieldVector<R, m> primitive_variables(0.);
+                const auto t = param.get("_t").at(0);
+                // density
+                primitive_variables[0] = 1. + 0.5 * std::sin(2. * M_PI * (t - 0.25));
+                // velocity
+                for (size_t ii = 0; ii < d; ++ii)
+                  primitive_variables[1 + ii] = 0.25 + 0.25 * std::sin(2. * M_PI * (t - 0.25));
+                // pressure
+                primitive_variables[m - 1] = 0.4;
+                return euler_tools.to_conservative(primitive_variables);
+              },
+              /*name=*/"periodic_density_variation",
+              /*parameter_type=*/{"_t", 1}));
+      local_periodic_density_variation_ = periodic_density_variation_->template as_grid_function<E>().local_function();
+      // inflow/outflow left, impermeable wall right
+      boundary_info = std::make_unique<XT::Grid::NormalBasedBoundaryInfo<XT::Grid::extract_intersection_t<GV>>>();
+      boundary_info->register_new_normal({-1}, new XT::Grid::InflowOutflowBoundary());
+      boundary_info->register_new_normal({1}, new XT::Grid::ImpermeableBoundary());
+      const XT::Grid::ApplyOn::CustomBoundaryIntersections<GV> impermeable_wall_filter(
+          *boundary_info, new XT::Grid::ImpermeableBoundary());
+      const XT::Grid::ApplyOn::CustomBoundaryIntersections<GV> inflow_outflow_filter(
+          *boundary_info, new XT::Grid::InflowOutflowBoundary());
+      auto op = std::make_unique<AdvectionFvOperator<GV, V, m>>(
+          space.grid_view(),
+          *self.numerical_flux_,
+          space,
+          space,
+          /*periodicity_restriction=*/*(inflow_outflow_filter || impermeable_wall_filter));
+      // the actual handling of inflow/outflow, see [DF2015, p. 421, (8.88)]
+      // (supposedly this does not work well for slow flows!)
+      const auto heuristic_euler_inflow_outflow_treatment = [&](const auto& intersection,
+                                                                const auto& xx_in_reference_intersection_coordinates,
+                                                                const auto& /*flux*/,
+                                                                const auto& u,
+                                                                const auto& param) {
+        using RangeType = XT::Common::FieldVector<R, m>;
+        // evaluate boundary values
+        const auto element = intersection.inside();
+        const auto xx_in_reference_element_coordinates =
+            intersection.geometryInInside().global(xx_in_reference_intersection_coordinates);
+        local_periodic_density_variation_->bind(element);
+        const RangeType bv = local_periodic_density_variation_->evaluate(xx_in_reference_element_coordinates, param);
+        // determine flow regime
+        const auto a = euler_tools.speed_of_sound_from_conservative(u);
+        const auto velocity = euler_tools.velocity_from_conservative(u);
+        const auto normal = intersection.unitOuterNormal(xx_in_reference_intersection_coordinates);
+        const auto flow_speed = velocity * normal;
+        // compute v
+        if (flow_speed < -a) {
+          // supersonic inlet
+          return bv;
+        } else if (!(flow_speed > 0)) {
+          // subsonic inlet
+          const auto rho_outer = euler_tools.density_from_conservative(bv);
+          const auto v_outer = euler_tools.velocity_from_conservative(bv);
+          const auto p_inner = euler_tools.pressure_from_conservative(u);
+          return euler_tools.to_conservative(XT::Common::hstack(rho_outer, v_outer, p_inner));
+        } else if (flow_speed < a) {
+          // subsonic outlet
+          const auto rho_inner = euler_tools.density_from_conservative(u);
+          const auto v_inner = euler_tools.velocity_from_conservative(u);
+          const auto p_outer = euler_tools.pressure_from_conservative(bv);
+          return euler_tools.to_conservative(XT::Common::hstack(rho_inner, v_inner, p_outer));
+        } else {
+          // supersonic outlet
+          return RangeType(u);
+        }
+      }; // ... heuristic_euler_inflow_outflow_treatment(...)
+      op->append(heuristic_euler_inflow_outflow_treatment, {}, inflow_outflow_filter);
+      // the actual handling of impermeable walls, see above
+      op->append(
+          /*boundary_extrapolation=*/
+          [&](const auto& intersection,
+              const auto& xx_in_reference_intersection_coordinates,
+              const auto& /*flux*/,
+              const auto& u,
+              const auto& /*param*/) {
             const auto normal = intersection.unitOuterNormal(xx_in_reference_intersection_coordinates);
             const auto rho = euler_tools.density_from_conservative(u);
             auto velocity = euler_tools.velocity_from_conservative(u);
@@ -203,6 +327,8 @@ protected:
   size_t visualization_steps_;
   std::string boundary_treatment;
   std::unique_ptr<XT::Grid::NormalBasedBoundaryInfo<XT::Grid::extract_intersection_t<GV>>> boundary_info;
+  std::unique_ptr<XT::Functions::LambdaFunction<d, m>> periodic_density_variation_;
+  std::unique_ptr<typename XT::Functions::ElementFunctionInterface<E, m>> local_periodic_density_variation_;
 }; // class InviscidCompressibleFlowEulerExplicitFvTest
 
 
@@ -226,12 +352,11 @@ TEST_F(InviscidCompressibleFlow1dEulerExplicitFvTest, impermeable_walls_by_invis
   this->boundary_treatment = "impermeable_walls_by_inviscid_mirror_treatment";
   this->run();
 }
-
-
-#if 0
-// TEST_F(InviscidCompressibleFlowEuler1dExplicitFV,
-//       inflow_from_the_left_by_heuristic_euler_treatment_impermeable_wall_right)
-//{
-//  this->inflow_from_the_left_by_heuristic_euler_treatment_impermeable_wall_right(10*100);
-//}
-#endif // 0
+TEST_F(InviscidCompressibleFlow1dEulerExplicitFvTest,
+       inflow_from_the_left_by_heuristic_euler_treatment_impermeable_wall_right)
+{
+  this->T_end_ = 2; // We need more time to hit the right wall
+  this->set_numerical_flux("vijayasundaram");
+  this->boundary_treatment = "inflow_from_the_left_by_heuristic_euler_treatment_impermeable_wall_right";
+  this->run();
+}