Add Functions::RayleighKotheVortex

doc/news/changes/minor/20230820Munch

@@ -0,0 +1,5 @@
+New: The Rayleigh--Kothe vortex has been extracted
+from step-68 and is now availabe as the new class
+Functions::RayleighKotheVortex.
+<br>
+(Bruno Blais, Peter Munch, 2023/08/20)

examples/step-68/step-68.cc

@@ -22,6 +22,7 @@
 #include <deal.II/base/bounding_box.h>
 #include <deal.II/base/conditional_ostream.h>
 #include <deal.II/base/discrete_time.h>
+#include <deal.II/base/function_lib.h>
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/parameter_acceptor.h>
 #include <deal.II/base/timer.h>
@@ -166,49 +167,6 @@ namespace Step68
 
 
 
-  // @sect3{Velocity profile}
-
-  // The velocity profile is provided as a Function object.
-  // This function is hard-coded within
-  // the example.
-  template <int dim>
-  class Vortex : public Function<dim>
-  {
-  public:
-    Vortex()
-      : Function<dim>(dim)
-    {}
-
-
-    virtual void vector_value(const Point<dim> &point,
-                              Vector<double>   &values) const override;
-  };
-
-
-  // The velocity profile for the Rayleigh-Kothe vertex is time-dependent.
-  // Consequently, the current time in the
-  // simulation (t) must be gathered from the Function object.
-  template <int dim>
-  void Vortex<dim>::vector_value(const Point<dim> &point,
-                                 Vector<double>   &values) const
-  {
-    const double T = 4;
-    const double t = this->get_time();
-
-    const double px = numbers::PI * point(0);
-    const double py = numbers::PI * point(1);
-    const double pt = numbers::PI / T * t;
-
-    values[0] = -2 * cos(pt) * pow(sin(px), 2) * sin(py) * cos(py);
-    values[1] = 2 * cos(pt) * pow(sin(py), 2) * sin(px) * cos(px);
-    if (dim == 3)
-      {
-        values[2] = 0;
-      }
-  }
-
-
-
   // @sect3{The <code>ParticleTracking</code> class declaration}
 
   // We are now ready to introduce the main class of our tutorial program.
@@ -278,7 +236,7 @@ namespace Step68
     MappingQ1<dim>                             mapping;
     LinearAlgebra::distributed::Vector<double> velocity_field;
 
-    Vortex<dim> velocity;
+    Functions::RayleighKotheVortex<dim> velocity;
 
     ConditionalOStream pcout;
 
@@ -305,6 +263,7 @@ namespace Step68
     , background_triangulation(mpi_communicator)
     , fluid_dh(background_triangulation)
     , fluid_fe(FE_Q<dim>(par.velocity_degree) ^ dim)
+    , velocity(4.0)
     , pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
     , interpolated_velocity(interpolated_velocity)
 

include/deal.II/base/function_lib.h

@@ -1752,6 +1752,72 @@ namespace Functions
     const std::vector<double> coefficients;
   };
 
+
+  /**
+   * A class that represents a time-dependent function object for a
+   * Rayleigh--Kothe vortex vector field. This is generally used as
+   * flow pattern in complex test cases for interface tracking methods
+   * (e.g., volume-of-fluid and level-set approaches) since it leads
+   * to strong rotation and elongation of the fluid @cite Blais2013.
+   *
+   * The stream function $\Psi$ of this Rayleigh-Kothe vortex is defined as:
+@f[
+\Psi = \frac{1}{\pi} \sin^2 (\pi x) \sin^2 (\pi y) \cos \left( \pi \frac{t}{T}
+\right)
+@f]
+   * where $T$ is half the period of the flow. The velocity profile in 2D
+($\textbf{u}=[u,v]^T$) is :
+@f{eqnarray*}{
+   u &=&  - \frac{\partial\Psi}{\partial y} = -2 \sin^2 (\pi x) \sin (\pi y)
+\cos (\pi y)  \cos \left( \pi \frac{t}{T} \right)\\ v &=&
+\frac{\partial\Psi}{\partial x} = 2 \cos(\pi x) \sin(\pi x) \sin^2 (\pi y) \cos
+\left( \pi \frac{t}{T} \right)
+@f}
+   * where $T$ is half the period of the flow.
+   *
+   * The velocity profile is illustrated in the following animation:
+   *
+@htmlonly
+<p align="center">
+  <iframe width="560" height="500"
+src="https://www.youtube.com/embed/m6hQm7etji8" frameborder="0"
+   allow="accelerometer; autoplay; encrypted-media; gyroscope;
+picture-in-picture" allowfullscreen></iframe>
+ </p>
+@endhtmlonly
+   *
+   * It can be seen that this velocity reverses periodically due to the term
+   * $\cos \left( \pi \frac{t}{T} \right)$ and that material will end up at its
+   * starting position after every period of length $t=2T$.
+   *
+   * @note This class is only implemented for 2D and 3D. In 3D, the
+   * third component is set to zero.
+   *
+   * @ingroup functions
+   */
+  template <int dim>
+  class RayleighKotheVortex : public Function<dim>
+  {
+  public:
+    /**
+     * Constructor.
+     */
+    RayleighKotheVortex(const double T = 1.0);
+
+    /**
+     * Return all components of a vector-valued function at a given point.
+     */
+    virtual void
+    vector_value(const Point<dim> &point,
+                 Vector<double>   &values) const override;
+
+  private:
+    /**
+     * Half the period of the flow.
+     */
+    const double T;
+  };
+
 #ifndef DOXYGEN
 
 

source/base/function_lib.cc

@@ -2947,6 +2947,32 @@ namespace Functions
            sizeof(coefficients);
   }
 
+  template <int dim>
+  RayleighKotheVortex<dim>::RayleighKotheVortex(const double T)
+    : Function<dim>(dim)
+    , T(T)
+  {
+    AssertThrow(dim > 1, ExcNotImplemented());
+  }
+
+
+  template <int dim>
+  void
+  RayleighKotheVortex<dim>::vector_value(const Point<dim> &point,
+                                         Vector<double>   &values) const
+  {
+    const double pi_x = numbers::PI * point(0);
+    const double pi_y = numbers::PI * point(1);
+    const double pi_t = numbers::PI / T * this->get_time();
+
+    values[0] = -2 * std::cos(pi_t) * std::pow(std::sin(pi_x), 2) *
+                std::sin(pi_y) * std::cos(pi_y);
+    values[1] = +2 * std::cos(pi_t) * std::pow(std::sin(pi_y), 2) *
+                std::sin(pi_x) * std::cos(pi_x);
+
+    if (dim == 3)
+      values[2] = 0;
+  }
 
 
   // explicit instantiations
@@ -3003,6 +3029,9 @@ namespace Functions
   template class Polynomial<1>;
   template class Polynomial<2>;
   template class Polynomial<3>;
+  template class RayleighKotheVortex<1>;
+  template class RayleighKotheVortex<2>;
+  template class RayleighKotheVortex<3>;
 } // namespace Functions
 
 DEAL_II_NAMESPACE_CLOSE

tests/performance/timing_step_68.cc

@@ -27,6 +27,7 @@
 #include <deal.II/base/bounding_box.h>
 #include <deal.II/base/conditional_ostream.h>
 #include <deal.II/base/discrete_time.h>
+#include <deal.II/base/function_lib.h>
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/parameter_acceptor.h>
 #include <deal.II/base/timer.h>
@@ -95,42 +96,6 @@ namespace Step68
   };
 
 
-  template <int dim>
-  class Vortex : public Function<dim>
-  {
-  public:
-    Vortex()
-      : Function<dim>(dim)
-    {}
-
-
-    virtual void
-    vector_value(const Point<dim> &point,
-                 Vector<double>   &values) const override;
-  };
-
-
-  template <int dim>
-  void
-  Vortex<dim>::vector_value(const Point<dim> &point,
-                            Vector<double>   &values) const
-  {
-    const double T = 4;
-    const double t = this->get_time();
-
-    const double px = numbers::PI * point(0);
-    const double py = numbers::PI * point(1);
-    const double pt = numbers::PI / T * t;
-
-    values[0] = -2 * cos(pt) * pow(sin(px), 2) * sin(py) * cos(py);
-    values[1] = 2 * cos(pt) * pow(sin(py), 2) * sin(px) * cos(px);
-    if (dim == 3)
-      {
-        values[2] = 0;
-      }
-  }
-
-
 
   template <int dim>
   class ParticleTracking
@@ -179,7 +144,7 @@ namespace Step68
     MappingQ1<dim>                             mapping;
     LinearAlgebra::distributed::Vector<double> velocity_field;
 
-    Vortex<dim> velocity;
+    Functions::RayleighKotheVortex<dim> velocity;
 
     ConditionalOStream pcout;
 
@@ -196,6 +161,7 @@ namespace Step68
     , background_triangulation(mpi_communicator)
     , fluid_dh(background_triangulation)
     , fluid_fe(FE_Q<dim>(par.velocity_degree), dim)
+    , velocity(4.0)
     , pcout(std::cout,
             debug && Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
     , interpolated_velocity(interpolated_velocity)