NonlinearLensInvariants.H

Go to the documentation of this file.

/* Copyright 2022-2023 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 IMPACTX_INVARIANTS_H
#define IMPACTX_INVARIANTS_H
 
#include "particles/ImpactXParticleContainer.H"
 
#include <AMReX_Extension.H>
#include <AMReX_REAL.H>
#include <AMReX_GpuComplex.H>
#include <AMReX_Math.H>
 
#include <cmath>
 
 
namespace impactx::diagnostics
{
    struct NonlinearLensInvariants
    {
        using PType = ImpactXParticleContainer::ParticleType;

        struct Data {
            amrex::ParticleReal H; 
            amrex::ParticleReal I; 
        };

        NonlinearLensInvariants (
            amrex::ParticleReal const alpha,
            amrex::ParticleReal const beta,
            amrex::ParticleReal const tn,
            amrex::ParticleReal const cn )
        : m_alpha(alpha), m_beta(beta), m_tn(tn), m_cn(cn)
        {
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        Data operator() (
            amrex::ParticleReal const x,
            amrex::ParticleReal const y,
            amrex::ParticleReal const px,
            amrex::ParticleReal const py
        ) const
        {
            using namespace amrex::literals; // for _rt and _prt
            using amrex::Math::powi;
 
            // a complex type with two amrex::ParticleReal
            using Complex = amrex::GpuComplex<amrex::ParticleReal>;
 
            // convert transverse phase space coordinates to normalized units
            amrex::ParticleReal const xn = x/(m_cn*std::sqrt(m_beta));
            amrex::ParticleReal const yn = y/(m_cn*std::sqrt(m_beta));
            amrex::ParticleReal const pxn = px*std::sqrt(m_beta)/m_cn + m_alpha*xn;
            amrex::ParticleReal const pyn = py*std::sqrt(m_beta)/m_cn + m_alpha*yn;
 
            // assign complex position zeta = x + iy
            Complex const zeta(xn, yn);
            Complex const zetaconj(xn, -yn);
            Complex const re1(1.0_prt, 0.0_prt);
            Complex const im1(0.0_prt, 1.0_prt);
 
            // compute croot = sqrt(1-zeta**2)
            Complex croot = powi<2>(zeta);
            croot = re1 - croot;
            croot = amrex::sqrt(croot);
 
            // compute carcsin = arcsin(zeta)
            Complex carcsin = im1*zeta + croot;
            carcsin = -im1*amrex::log(carcsin);
 
            // compute complex potential appearing in H invariant
            Complex Hpotential = zeta/croot;
            Hpotential = Hpotential*carcsin;
 
            // compute complex potential appearing in I invariant
            Complex Ipotential = (zeta+zetaconj)/croot;
            Ipotential = Ipotential*carcsin;
 
            // evaluate real parts
            amrex::ParticleReal Hinv = Hpotential.m_real;
            amrex::ParticleReal Iinv = Ipotential.m_real;
 
            // compute invariants H and I
            amrex::ParticleReal const Jz = xn*pyn - yn*pxn;
            Hinv = (powi<2>(xn) + powi<2>(yn) + powi<2>(pxn) + powi<2>(pyn))/2
                 + m_tn*Hinv;
            Iinv = powi<2>(Jz) + powi<2>(pxn) + powi<2>(xn) + m_tn*Iinv;
 
            return {Hinv, Iinv};
 
        }
 
    private:
        amrex::ParticleReal m_alpha; 
        amrex::ParticleReal m_beta; 
        amrex::ParticleReal m_tn; 
        amrex::ParticleReal m_cn; 
    };
 
} // namespace impactx
 
#endif // IMPACTX_INVARIANTS_H