From: Peter Munch Date: Sun, 20 Aug 2023 09:31:21 +0000 (+0200) Subject: Add Functions::RayleighKotheVortex X-Git-Tag: relicensing~531^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15932%2Fhead;p=dealii.git Add Functions::RayleighKotheVortex --- diff --git a/doc/news/changes/minor/20230820Munch b/doc/news/changes/minor/20230820Munch new file mode 100644 index 0000000000..67cb25ff9f --- /dev/null +++ b/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. +
+(Bruno Blais, Peter Munch, 2023/08/20) diff --git a/examples/step-68/step-68.cc b/examples/step-68/step-68.cc index a91c59e197..d47649f3bf 100644 --- a/examples/step-68/step-68.cc +++ b/examples/step-68/step-68.cc @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -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 - class Vortex : public Function - { - public: - Vortex() - : Function(dim) - {} - - - virtual void vector_value(const Point &point, - Vector &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 - void Vortex::vector_value(const Point &point, - Vector &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 ParticleTracking class declaration} // We are now ready to introduce the main class of our tutorial program. @@ -278,7 +236,7 @@ namespace Step68 MappingQ1 mapping; LinearAlgebra::distributed::Vector velocity_field; - Vortex velocity; + Functions::RayleighKotheVortex velocity; ConditionalOStream pcout; @@ -305,6 +263,7 @@ namespace Step68 , background_triangulation(mpi_communicator) , fluid_dh(background_triangulation) , fluid_fe(FE_Q(par.velocity_degree) ^ dim) + , velocity(4.0) , pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0) , interpolated_velocity(interpolated_velocity) diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h index 3a9baf46db..075c6541c7 100644 --- a/include/deal.II/base/function_lib.h +++ b/include/deal.II/base/function_lib.h @@ -1752,6 +1752,72 @@ namespace Functions const std::vector 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 +

+ +

+@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 + class RayleighKotheVortex : public Function + { + 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 &point, + Vector &values) const override; + + private: + /** + * Half the period of the flow. + */ + const double T; + }; + #ifndef DOXYGEN diff --git a/source/base/function_lib.cc b/source/base/function_lib.cc index 0209f542c9..221ceaaa35 100644 --- a/source/base/function_lib.cc +++ b/source/base/function_lib.cc @@ -2947,6 +2947,32 @@ namespace Functions sizeof(coefficients); } + template + RayleighKotheVortex::RayleighKotheVortex(const double T) + : Function(dim) + , T(T) + { + AssertThrow(dim > 1, ExcNotImplemented()); + } + + + template + void + RayleighKotheVortex::vector_value(const Point &point, + Vector &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 diff --git a/tests/performance/timing_step_68.cc b/tests/performance/timing_step_68.cc index fdff901c76..565580f60e 100644 --- a/tests/performance/timing_step_68.cc +++ b/tests/performance/timing_step_68.cc @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -95,42 +96,6 @@ namespace Step68 }; - template - class Vortex : public Function - { - public: - Vortex() - : Function(dim) - {} - - - virtual void - vector_value(const Point &point, - Vector &values) const override; - }; - - - template - void - Vortex::vector_value(const Point &point, - Vector &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 class ParticleTracking @@ -179,7 +144,7 @@ namespace Step68 MappingQ1 mapping; LinearAlgebra::distributed::Vector velocity_field; - Vortex velocity; + Functions::RayleighKotheVortex velocity; ConditionalOStream pcout; @@ -196,6 +161,7 @@ namespace Step68 , background_triangulation(mpi_communicator) , fluid_dh(background_triangulation) , fluid_fe(FE_Q(par.velocity_degree), dim) + , velocity(4.0) , pcout(std::cout, debug && Utilities::MPI::this_mpi_process(mpi_communicator) == 0) , interpolated_velocity(interpolated_velocity)