ERF_MorrisonGammaFunction.H File Reference

#include <string>
#include <vector>
#include <memory>
#include <complex>
#include <cmath>
#include <AMReX_Math.H>
#include <AMReX_FArrayBox.H>
#include <AMReX_Geometry.H>
#include <AMReX_TableData.H>
#include <AMReX_MultiFabUtil.H>
#include "ERF.H"
#include "ERF_Constants.H"
#include "ERF_MicrophysicsUtils.H"
#include "ERF_IndexDefines.H"
#include "ERF_DataStruct.H"
#include "ERF_NullMoist.H"
#include "ERF_Morrison.H"

Include dependency graph for ERF_MorrisonGammaFunction.H:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	wrf_gamma (amrex::Real x)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	gamma_function (Real x)

Variables
constexpr Real	xxx = Real(0.9189385332046727417803297)

Function Documentation

gamma_function()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real gamma_function

(

Real

x

)

                            {
  return wrf_gamma(x);
}

Referenced by Morrison::Advance().

Here is the call graph for this function:

Here is the caller graph for this function:

wrf_gamma()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real wrf_gamma

(

amrex::Real

x

)

{
    // Debug: using printf since it's GPU compatible
//    printf("wrf_gamma: Input value x = %g\n", x);
 
    // Local variables
    int i, n;
    bool parity = false;
    amrex::Real fact, res, sum, twelve, xbig, xden, xinf, xminin;
    amrex::Real xnum, y, y1, ysq, z;
    amrex::Real c[7];
    amrex::Real p[8];
    amrex::Real q[8];
 
    // Mathematical constants
    twelve = Real(12.0);
 
    // Machine dependent parameters
    xbig = Real(35.040);
    xminin = Real(1.18e-38);
    amrex::Real eps = Real(1.19e-7);
    xinf = Real(3.4e38);
 
    // Numerator and denominator coefficients for rational minimax approximation over (1,2)
    p[0] = -Real(1.71618513886549492533811e+0);
    p[1] =  Real(2.47656508055759199108314e+1);
    p[2] = -Real(3.79804256470945635097577e+2);
    p[3] =  Real(6.29331155312818442661052e+2);
    p[4] =  Real(8.66966202790413211295064e+2);
    p[5] = -Real(3.14512729688483675254357e+4);
    p[6] = -Real(3.61444134186911729807069e+4);
    p[7] =  Real(6.64561438202405440627855e+4);
 
    q[0] = -Real(3.08402300119738975254353e+1);
    q[1] =  Real(3.15350626979604161529144e+2);
    q[2] = -Real(1.01515636749021914166146e+3);
    q[3] = -Real(3.10777167157231109440444e+3);
    q[4] =  Real(2.25381184209801510330112e+4);
    q[5] =  Real(4.75584627752788110767815e+3);
    q[6] = -Real(1.34659959864969306392456e+5);
    q[7] = -Real(1.15132259675553483497211e+5);
 
    // Coefficients for minimax approximation over (12, inf)
    c[0] = -Real(1.910444077728e-03);
    c[1] = Real(8.4171387781295e-04);
    c[2] = -Real(5.952379913043012e-04);
    c[3] = Real(7.93650793500350248e-04);
    c[4] = -Real(2.777777777777681622553e-03);
    c[5] = Real(8.333333333333333331554247e-02);
    c[6] = Real(5.7083835261e-03);
 
    // Initialize variables
    parity = false;
    fact = one;
    n = 0;
    y = x;
 
//    printf("wrf_gamma: Initial y = %g\n", y);
 
    if (y <= Real(0)) {
        // Argument is negative
//        printf("wrf_gamma: Handling negative argument\n");
        y = -x;
        y1 = std::floor(y);
        res = y - y1;
        if (res != Real(0)) {
            if (y1 != std::floor(y1 * myhalf) * Real(2))
                parity = true;
            Real pi=amrex::Math::pi<Real>();
            fact = -pi / std::sin(pi * res);
            y = y + one;
//            printf("wrf_gamma: After reflection formula: y = %g, fact = %g, parity = %d\n",
//                   y, fact, parity);
        }
        else {
//            printf("wrf_gamma: Singularity detected, returning xinf = %g\n", xinf);
            res = xinf;
            return res;
        }
    }
 
    // Argument is positive
    if (y < eps) {
        // Argument < eps
//        printf("wrf_gamma: Small argument branch (y < eps)\n");
        if (y >= xminin) {
            res = one / y;
//            printf("wrf_gamma: Small argument result: res = %g\n", res);
        }
        else {
//            printf("wrf_gamma: Argument too small, returning xinf = %g\n", xinf);
            res = xinf;
            return res;
        }
    }
    else if (y < twelve) {
        // Medium range argument
//        printf("wrf_gamma: Medium range branch (eps <= y < 12)\n");
        y1 = y;
        if (y < one) {
            // Real(0) < argument < one
//            printf("wrf_gamma: Sub-branch: 0 < y < 1\n");
            z = y;
            y = y + one;
        }
        else {
            // one < argument < Real(12.0), reduce argument if necessary
            n = static_cast<int>(y) - 1;
//            printf("wrf_gamma: Sub-branch: 1 <= y < 12, n = %d\n", n);
            y = y - static_cast<amrex::Real>(n);
            z = y - one;
        }
 
        // Evaluate approximation
//        printf("wrf_gamma: Before approximation: z = %g, y = %g\n", z, y);
        xnum = Real(0);
        xden = one;
        for (i = 0; i < 8; i++) {
            xnum = (xnum + p[i]) * z;
            xden = xden * z + q[i];
        }
        res = xnum / xden + one;
//        printf("wrf_gamma: After approximation: res = %g\n", res);
 
        if (y1 < y) {
            // Adjust result for case Real(0) < argument < one
            res = res / y1;
//            printf("wrf_gamma: Adjusted for y < 1: res = %g\n", res);
        }
        else if (y1 > y) {
            // Adjust for two < argument < Real(12.0)
//            printf("wrf_gamma: Adjusting for y > 2 with %d multiplications\n", n);
            for (i = 0; i < n; i++) {
                res = res * y;
                y = y + one;
//                printf("wrf_gamma: Multiplication %d: res = %g, y = %g\n", i+1, res, y);
            }
        }
    }
    else {
        // Large argument
//        printf("wrf_gamma: Large argument branch (y >= 12)\n");
        if (y <= xbig) {
            ysq = y * y;
            sum = c[6];
            for (i = 0; i < 6; i++) {
                sum = sum / ysq + c[i];
//                printf("wrf_gamma: Sum step %d: sum = %g\n", i+1, sum);
            }
            sum = sum / y - y + xxx;
            sum = sum + (y - myhalf) * std::log(y);
//            printf("wrf_gamma: Before exp: sum = %g\n", sum);
            res = std::exp(sum);
//            printf("wrf_gamma: After exp: res = %g\n", res);
        }
        else {
//            printf("wrf_gamma: Argument too large, returning xinf = %g\n", xinf);
            res = xinf;
            return res;
        }
    }
 
    // Final adjustments
    if (parity) {
        res = -res;
//        printf("wrf_gamma: Applied parity adjustment: res = %g\n", res);
    }
    if (fact != one) {
        res = fact / res;
//        printf("wrf_gamma: Applied reflection adjustment: res = %g\n", res);
    }
 
//    printf("wrf_gamma: Final result = %g\n", res);
    return res;
}

Referenced by gamma_function().

Here is the caller graph for this function:

Variable Documentation

xxx

constexpr Real xxx = Real(0.9189385332046727417803297)

constexpr

Referenced by wrf_gamma().