Implementation of Bessel functions of second kind (yn)

Docs/sphinx_documentation/source/Basics.rst

@@ -595,7 +595,7 @@ numbers can be computed with ``min`` and ``max``, respectively.  It supports
 the Heaviside step function, ``heaviside(x1,x2)`` that gives ``0``, ``x2``,
 ``1``, for ``x1 < 0``, ``x1 = 0`` and ``x1 > 0``, respectively.
 It supports the Bessel function of the first kind of order ``n``
-``jn(n,x)``. Complete elliptic integrals of the first and second kind, ``comp_ellint_1(k)`` and ``comp_ellint_2(k)``,
+``jn(n,x)``, and the Bessel function of the second kind of order ``n`` ``yn(n,x)``. Complete elliptic integrals of the first and second kind, ``comp_ellint_1(k)`` and ``comp_ellint_2(k)``,
 are supported.
 There is ``if(a,b,c)`` that gives ``b`` or ``c`` depending on the value of
 ``a``.  A number of comparison operators are supported, including ``<``,

Src/Base/Parser/AMReX_Parser_Y.H

@@ -91,6 +91,7 @@ enum parser_f2_t {  // Built-in functions with two arguments
     PARSER_OR,
     PARSER_HEAVISIDE,
     PARSER_JN,
+    PARSER_YN,
     PARSER_MIN,
     PARSER_MAX,
     PARSER_FMOD
@@ -116,6 +117,7 @@ std::string_view parser_f2_s[] =
     "or",
     "heaviside",
     "jn",
+    "yn",
     "min",
     "max",
     "fmod"
@@ -400,6 +402,21 @@ T parser_math_jn (int a, T b)
 #endif
 }
 
+template <typename T>
+AMREX_GPU_HOST_DEVICE AMREX_NO_INLINE
+T parser_math_yn (int a, T b)
+{
+#if defined AMREX_USE_SYCL || defined __MINGW32__
+    amrex::ignore_unused(a,b);
+    // neither yn(f) nor std::cyl_bessel_y work yet
+    // https://github.com/oneapi-src/oneAPI-spec/issues/308
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "parser: yn in SYCL not supported yet");
+    return 0.0;
+#else
+    return yn(a, b);
+#endif
+}
+
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE double
 parser_call_f1 (enum parser_f1_t type, double a)
 {
@@ -460,6 +477,8 @@ parser_call_f2 (enum parser_f2_t type, double a, double b)
         return (a < 0.0) ? 0.0 : ((a > 0.0) ? 1.0 : b);
     case PARSER_JN:
         return parser_math_jn<double>(int(a),b);
+    case PARSER_YN:
+        return parser_math_yn<double>(int(a),b);
     case PARSER_MIN:
         return (a < b) ? a : b;
     case PARSER_MAX:

Src/Base/Parser/amrex_parser.l

@@ -70,6 +70,7 @@ EXP	([Ee][-+]?[0-9]+)
 "pow"         { amrex_parserlval.f2 = amrex::PARSER_POW;       return F2; }
 "heaviside"   { amrex_parserlval.f2 = amrex::PARSER_HEAVISIDE; return F2; }
 "jn"          { amrex_parserlval.f2 = amrex::PARSER_JN;        return F2; }
+"yn"          { amrex_parserlval.f2 = amrex::PARSER_YN;        return F2; }
 "min"         { amrex_parserlval.f2 = amrex::PARSER_MIN;       return F2; }
 "max"         { amrex_parserlval.f2 = amrex::PARSER_MAX;       return F2; }
 "fmod"        { amrex_parserlval.f2 = amrex::PARSER_FMOD;      return F2; }