CovarianceMatrixMath.H

Go to the documentation of this file.

/* Copyright 2022-2026 The Regents of the University of California, through Lawrence
 *           Berkeley National Laboratory (subject to receipt of any required
 *           approvals from the U.S. Dept. of Energy). All rights reserved.
 *
 * This file is part of ImpactX.
 *
 * Authors: Chad Mitchell, Axel Huebl
 * License: BSD-3-Clause-LBNL
 */
#ifndef COVARIANCE_MATRIX_MATH_H
#define COVARIANCE_MATRIX_MATH_H
 
#include <ablastr/constant.H>
#include <ablastr/warn_manager/WarnManager.H>
 
#include <AMReX_REAL.H>
#include <AMReX_GpuComplex.H>
#include <AMReX_Math.H>
 
#include <cmath>
#include <tuple>
 
 
namespace impactx::diagnostics
{

    std::tuple<
        amrex::ParticleReal,
        amrex::ParticleReal,
        amrex::ParticleReal
    >
    CubicRootsTrig (
        amrex::ParticleReal a,
        amrex::ParticleReal b,
        amrex::ParticleReal c,
        amrex::ParticleReal d
    )
    {
        using namespace amrex::literals;
        using amrex::Math::powi;
        using ablastr::constant::math::pi;
 
        std::tuple<amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal> roots;
        amrex::ParticleReal x1 = 0.0_prt;
        amrex::ParticleReal x2 = 0.0_prt;
        amrex::ParticleReal x3 = 0.0_prt;
 
        amrex::ParticleReal Q = (3.0_prt*a*c - powi<2>(b))/(9.0_prt * powi<2>(a));
        amrex::ParticleReal R = (9.0_prt*a*b*c - 27_prt*powi<2>(a)*d - 2.0_prt*powi<3>(b))/(54.0_prt*powi<3>(a));
        amrex::ParticleReal discriminant = powi<3>(Q) + powi<2>(R);
 
        // Discriminant should be < 0.  Otherwise, keep theta at default and throw an error.
        amrex::ParticleReal tol = 1.0e-12_prt;  //allow for roundoff error
        if (discriminant > tol) {
 
            ablastr::warn_manager::WMRecordWarning(
                    "Impactx::diagnostics::CubicRootsTrig",
                    "Polynomial appearing in CubicRootsTrig has one or more complex "
                    "(non-real) roots.  Only the real part is returned.  This "
                    "suggests a loss of numerical precision in computation of the "
                    "eigenemittances.  Treat eigenemittance values with caution.",
                    ablastr::warn_manager::WarnPriority::medium
            );
 
           std::cout << "Polynomial in CubicRoots has one or more complex roots." << "\n";
 
        } else if (Q == 0.0_prt) {  // Special case of a triple root
 
           x1 = - b/(3.0_prt*a);
           x2 = - b/(3.0_prt*a);
           x3 = - b/(3.0_prt*a);
 
        } else {
 
           //Three real roots in trigonometric form.
           amrex::ParticleReal theta = std::acos(R/std::sqrt(-powi<3>(Q)));
           x1 = 2.0_prt*std::sqrt(-Q)*std::cos(theta/3.0_prt) - b/(3.0_prt*a);
           x2 = 2.0_prt*std::sqrt(-Q)*std::cos(theta/3.0_prt + 2.0_prt*pi/3.0_prt) - b/(3.0_prt*a);
           x3 = 2.0_prt*std::sqrt(-Q)*std::cos(theta/3.0_prt + 4.0_prt*pi/3.0_prt) - b/(3.0_prt*a);
 
        }
 
        roots = std::make_tuple(x1,x2,x3);
        return roots;
    }
 

    std::tuple<
        amrex::ParticleReal,
        amrex::ParticleReal,
        amrex::ParticleReal
    >
    CubicRootsAlg (
        amrex::ParticleReal a,
        amrex::ParticleReal b,
        amrex::ParticleReal c,
        amrex::ParticleReal d
    )
    {
        using namespace amrex::literals;
        using amrex::Math::powi;
        using Complex = amrex::GpuComplex<amrex::ParticleReal>;
 
        std::tuple<amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal> roots;
        amrex::ParticleReal x1 = 0_prt;
        amrex::ParticleReal x2 = 0_prt;
        amrex::ParticleReal x3 = 0_prt;
 
        amrex::ParticleReal Q = (3_prt*a*c - powi<2>(b))/(9_prt * powi<2>(a));
        amrex::ParticleReal R = (9_prt*a*b*c - 27_prt*powi<2>(a)*d - 2_prt*powi<3>(b))/(54_prt*powi<3>(a));
        amrex::ParticleReal discriminant = powi<3>(Q) + powi<2>(R);
 
        // Define complex variable C
        Complex Qc(Q,0_prt);
        Complex Rc(R,0_prt);
        Complex Dc(discriminant,0_prt);
        Complex C = amrex::pow(-Rc + amrex::sqrt(Dc),1_prt/3_prt);
 
        // Define a cubic root of unity xi
        amrex::ParticleReal xire = -0.5_prt;
        amrex::ParticleReal xiim = std::sqrt(3_prt)*0.5_prt;
        Complex xi(xire,xiim);
 
        //Three roots in algebraic form.
 
        if (C.m_real == 0_prt && C.m_imag == 0_prt) {  // Special case of a triple root
 
           x1 = - b/(3_prt*a);
           x2 = - b/(3_prt*a);
           x3 = - b/(3_prt*a);
 
        } else {
 
           Complex z1 = Qc/C - C;
           Complex z2 = Qc/(xi*C) - xi*C;
           Complex z3 = Qc/(powi<2>(xi)*C) - powi<2>(xi)*C;
           x1 = z2.m_real - b/(3_prt*a);
           x2 = z1.m_real - b/(3_prt*a);
           x3 = z3.m_real - b/(3_prt*a);
        }
 
        roots = std::make_tuple(x1,x2,x3);
        return roots;
    }
 
} // namespace impactx::diagnostics
 
#endif // COVARIANCE_MATRIX_MATH_H