add_f0_analytic.hpp

Go to the documentation of this file.

#ifndef ADD_F0_ANALYTIC_HPP
#define ADD_F0_ANALYTIC_HPP

#include <Kokkos_Core.hpp>
#include <string>
#include "space_settings.hpp"
#include "sml.hpp"
#include "plasma.hpp"
#include "electric_field.hpp"
#include "species.hpp"
#include "vgrid_distribution.hpp"

void add_f0_analytic(const Simulation<DeviceType>& sml, ElectricField<DeviceType>& electric_field,
                     const DomainDecomposition<DeviceType>& pol_decomp,
                     Plasma& plasma, VGridDistribution<HostType>& f0_f, bool subtract = false){

    double sign = (subtract ? -1.0 : 1.0);

    int nvr = f0_f.n_vr(), nvz = f0_f.n_vz();
    int nvp = (nvz - 1)/2;
    double dsmu = f0_f.dsmu, dvp = f0_f.dvp;
    double inv_mu0_factor = f0_f.inv_mu0_factor;
    plasma.for_all_nonadiabatic_species([&](Species<DeviceType>& species){
        int isp = species.nonadiabatic_idx;
        Kokkos::parallel_for("add_f0_analytic", Kokkos::RangePolicy<HostExSpace>( 0, f0_f.n_nodes()), KOKKOS_LAMBDA(const int node){
            // Add contribution from anomalous diffusion
            int node_global = node+pol_decomp.node_offset;
            double prefac2 = species.f0.temp_ev_h(node)/(species.f0.temp_ev_h(node) + species.f0.delta_T_h(node_global));
            double prefac1 = sqrt(prefac2)*(species.f0.den_h(node) + species.f0.delta_n_h(node_global))
                /(species.f0.temp_ev_h(node) + species.f0.delta_T_h(node_global));
            double flow = species.f0.flow_h(node) + species.f0.delta_u_h(node_global)/species.get_f0_eq_thermal_velocity_lnode_h(node);

            double pot = 0.0;
            if(species.is_electron){
                // normalized potential

                // Get dpot
                Simd<double> dpot;
                const int i_simd = 0;
                dpot.zero();
#ifdef XGCA
#ifdef COL_RELAX_TEST
                // Keep dpot[i_simd] = 0
#else
                // Upper limit for pot (We use linearized Poisson equation!!!)
                // This is important at the wall strike points
                electric_field.dpot_h(0, node_global).gather(dpot, i_simd, 1.0);
#endif
#elif defined(XGC1)
                // Old Fortran remnant

                // Average the two planes
                double wphi[2];
                wphi[0] = 0.5;
                wphi[1] = 0.5;
                electric_field.dpot_ff_h(node_global).gather(dpot, i_simd, 1.0, wphi);

                // Old Fortran remnant
#endif
                pot = dpot[i_simd]/(species.f0.temp_ev_h(node) + species.f0.delta_T_h(node_global));
                pot = max(-sml.dpot_te_limit, min(sml.dpot_te_limit, pot));
            }

            for (int ivr = 0; ivr < nvr; ivr++){
                // normalized
                double smu = ivr*dsmu; // mu or sqrt mu (normalized)
                double mu = smu*smu;   // now mu is normalized magnetic moment or (v_perp/v_th)^2

                for(int ivz = 0; ivz < nvz; ivz++){
                    // Shift ivz so it runs between -nvp and nvp
                    double vp = (ivz - nvp)*dvp - flow;

                    if(ivr == 0){
                        smu = dsmu*inv_mu0_factor;
                    }
                    double en = 0.5*(vp*vp + mu);
                    // 'smu' is v_perp/v_th

                    f0_f(isp, ivr, node, ivz) += sign*prefac1*smu*exp(-prefac2*en + pot);
                }
            }
        });
        Kokkos::fence();
    }, Plasma::NoDevicePtl);
}


void f0_remove_negative(Plasma& plasma, VGridDistribution<HostType>& f0_f){

    constexpr double eps=1e-2;

    int nvr = f0_f.n_vr(), nvz = f0_f.n_vz();
    int nvp = (nvz - 1)/2;
    double dsmu = f0_f.dsmu, dvp = f0_f.dvp;
    double inv_mu0_factor = f0_f.inv_mu0_factor;

    plasma.for_all_nonadiabatic_species([&](Species<DeviceType>& species){
        int isp = species.nonadiabatic_idx;
        Kokkos::parallel_for("f0_remove_negative", Kokkos::RangePolicy<HostExSpace>( 0, f0_f.n_nodes()), KOKKOS_LAMBDA(const int node){
            double den_temp = species.f0.n_Ta_h(node);
            double flow = species.f0.flow_h(node);
            for (int ivr = 0; ivr < nvr; ivr++){
                double smu = ivr*dsmu;
                double mu = smu*smu;
                if(ivr == 0){
                    smu = dsmu*inv_mu0_factor;
                }
                for(int ivz = 0; ivz < nvz; ivz++){
                    if(f0_f(isp, ivr, node, ivz) <= 0.0){
                        // Use <= to make sure we never have density = 0
                        // For the full distribution function, this
                        // is a small number
                        // Shift ivz so it runs between -nvp and nvp
                        double en = 0.5*( ((ivz - nvp)*dvp - flow)*
                            ((ivz - nvp)*dvp - flow) + mu);
                        f0_f(isp, ivr, node, ivz) = eps*smu*den_temp*exp(-en);
                    }
                }
            }
        });
        Kokkos::fence();
    }, Plasma::NoDevicePtl);
}
#endif