varying_background.hpp

Go to the documentation of this file.

#ifndef VARYING_BACKGROUND_HPP
#define VARYING_BACKGROUND_HPP

#include "magnetic_field.hpp"
#include "plasma.hpp"
#include "get_monte_num.hpp"
#include "get_volume.hpp"

template<class Device>
class VaryingBackground{

    enum{
        Density1=0,
        Energy,
        VPara,
        Density2,
        VPara2,
        NTmp
    };

    enum{
        Density=0,
        Temp,
        Flow,
        N
    };

    double pin;
    double pout;
    int m;
    int mtheta;
    double inv_dp;
    double inv_dtheta;

    View<double**, CLayout, Device> vol; 
    View<double****, CLayout, Device> bg; 

    public:

    int period; 

    VaryingBackground(){}

    VaryingBackground(NLReader::NamelistReader& nlr, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, int n_nonadiabatic_species){
        nlr.use_namelist("col_param");
        period = nlr.get<int>("col_vb_period", 1);
        pin = nlr.get<double>("col_vb_pin", magnetic_field.inpsi);
        pout = nlr.get<double>("col_vb_pout", std::min(magnetic_field.outpsi, 1.03));
        m = nlr.get<int>("col_vb_m", 50);
        mtheta = nlr.get<int>("col_vb_mtheta", 8);

        // Normalize
        if(nlr.present("col_vb_pin")) pin *= magnetic_field.psi_norm();
        if(nlr.present("col_vb_pout")) pout *= magnetic_field.psi_norm();

        double dp = (pout - pin)/m;
        double dtheta = TWOPI/mtheta;
        inv_dp = 1.0/dp;
        inv_dtheta = 1.0/dtheta;

        // Allocate views
        bg = View<double****, CLayout, Device>(NoInit("vb_bg"), n_nonadiabatic_species, m+1, mtheta, N);
        vol = View<double**, CLayout, Device>("vb_vol", m+1, mtheta);

        int n_monte_carlo = get_monte_num(nlr, grid.nnode, grid.nplanes); // Note: n_ptl is needed for sml_monte_num
        bool use_nonrandom_sampling = false;
        GPTLstart("COL_VB_VOL");
        monte_carlo_col_vb_vol(pol_decomp, grid, magnetic_field, n_monte_carlo, use_nonrandom_sampling,            
                               pin, pout, m, mtheta, inv_dp, inv_dtheta, vol);
        GPTLstart("COL_VB_VOL");
    }

    KOKKOS_INLINE_FUNCTION bool is_in_range(const MagneticField<DeviceType>& magnetic_field, double r, double z, double psi) const {
        return pin<psi && psi<pout && z>magnetic_field.equil.xpt.z && magnetic_field.is_in_region_1_or_2(r,z,psi);
    }

    KOKKOS_INLINE_FUNCTION void interp_bg_profile(double psi, double theta, int isp, double& den, double& temp, double& up) const {
        double tmp = (psi-pin)*inv_dp;
        int j = int(tmp);
        j = min(m-1,max(0,j));
        int l = ((int)(floor(theta*inv_dtheta + 0.5)))%mtheta;
        double bb = tmp - j;
        double aa = 1.0 - bb;

        den = bg(isp,j,l,Density)*aa + bg(isp,j+1,l,Density)*bb;
        temp = bg(isp,j,l,Temp)*aa + bg(isp,j+1,l,Temp)*bb;

        j = floor((psi-pin)*inv_dp + 0.5);
        up = bg(isp,j,l,Flow);
    }

    // get ion/electron density and temperature for collision routine
    void update(const MagneticField<DeviceType>& magnetic_field, Species<DeviceType>& species){
        View<double***, CLayout, DeviceType> dum(NoInit("dum"),m+1,mtheta,N);
        Kokkos::deep_copy(bg, 1.0e-99); // to avoid divide by zero
        Kokkos::deep_copy(dum, 1.0e-99); // to avoid divide by zero

        auto vb = *this; // Use a local shallow copy so that the lambda function knows to copy it

        species.for_all_particles("col_snapshot", KOKKOS_LAMBDA( const int idx ){
            // Get index info
            AoSoAIndices<Device> inds(idx);

            // Copy from AoSoA to SoA
            SimdParticles part;
            AoSoA_2_simd(part, species.ptl(), inds.a, inds.s);

            Simd<double> psi;
            magnetic_field.get_psi(part.ph.v(),psi);

            Simd<double> Bmag;
            magnetic_field.bmag_interpol(part.ph.v(), Bmag);

            Simd<double> theta;
            magnetic_field.get_theta(part.ph.x(), theta);

            FORCE_SIMD
            for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
                if(part.is_active(i_simd)){
                    if(vb.pin<psi[i_simd] && psi[i_simd]<vb.pout && magnetic_field.is_in_region_1_or_2(part.ph.r[i_simd],part.ph.z[i_simd],psi[i_simd])){
                        //energy calculation
                        double weight = part.ct.w0[i_simd];
                        double upar = part.ph.rho[i_simd]*Bmag[i_simd]*species.c_m;
                        double energy = 0.5*species.mass*upar*upar + part.ct.mu[i_simd]*Bmag[i_simd];

                        double tmp = (psi[i_simd] - vb.pin)*vb.inv_dp;
                        int j = (int)(tmp);
                        j = min(vb.m-1,max(0,j));
                        int l = ((int)(floor(theta[i_simd]*vb.inv_dtheta + 0.5)))%vb.mtheta;
                        double bb = tmp-j;
                        double aa = 1.0-bb;

                        Kokkos::atomic_add(&(dum(j,l,Density1)), weight*aa);
                        Kokkos::atomic_add(&(dum(j,l,Energy)), weight*energy*aa);
                        Kokkos::atomic_add(&(dum(j,l,VPara)), weight*upar*aa);
                        Kokkos::atomic_add(&(dum(j+1,l,Density1)), weight*bb);
                        Kokkos::atomic_add(&(dum(j+1,l,Energy)), weight*energy*bb);
                        Kokkos::atomic_add(&(dum(j+1,l,VPara)), weight*upar*bb);


                        j = floor((psi[i_simd]-vb.pin)*vb.inv_dp + 0.5);
                        Kokkos::atomic_add(&(dum(j,l,VPara2)), weight);
                        Kokkos::atomic_add(&(dum(j,l,Density2)), weight*upar);
                    }
                }
            }
        });

#ifdef USE_MPI
        TIMER("COL_SNAP_RED",
           MPI_Allreduce(MPI_IN_PLACE, dum.data(), dum.size(), MPI_DOUBLE, MPI_SUM, SML_COMM_WORLD) );
#endif

        int isp = species.nonadiabatic_idx;
        int n = (m+1)*mtheta;
        Kokkos::parallel_for("snapshot_convert", Kokkos::RangePolicy<ExSpace>(0, n), KOKKOS_LAMBDA( const int idx){
            int i = idx%vb.mtheta; // Faster index
            int j = idx/vb.mtheta;
            constexpr double prefac = 2.0/3.0*J_2_EV; // Convert to eV
            const double vpara_normed = dum(j,i,VPara)/dum(j,i,Density1);
            const double epara = 0.5*species.mass*(vpara_normed*vpara_normed);
            vb.bg(isp,j,i,Temp) = prefac*(dum(j,i,Energy) - epara)/dum(j,i,Density1);
            vb.bg(isp,j,i,Density) = dum(j,i,Density1)/vb.vol(j,i);
            vb.bg(isp,j,i,Flow) = dum(j,i,VPara2)/dum(j,i,Density2); // mean flow for moving frame collision
        });
    }
};

#endif