XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
add_f0_analytic.hpp
Go to the documentation of this file.
1 #ifndef ADD_F0_ANALYTIC_HPP
2 #define ADD_F0_ANALYTIC_HPP
3 
4 #include <Kokkos_Core.hpp>
5 #include <string>
6 #include "space_settings.hpp"
7 #include "sml.hpp"
8 #include "plasma.hpp"
9 #include "electric_field.hpp"
10 #include "species.hpp"
11 #include "vgrid_distribution.hpp"
12 
26  const DomainDecomposition<DeviceType>& pol_decomp,
27  Plasma& plasma, VGridDistribution<HostType>& f0_f, bool subtract = false){
28 
29  double sign = (subtract ? -1.0 : 1.0);
30 
31  int nvr = f0_f.n_vr(), nvz = f0_f.n_vz();
32  int nvp = (nvz - 1)/2;
33  double dsmu = f0_f.dsmu, dvp = f0_f.dvp;
34  double inv_mu0_factor = f0_f.inv_mu0_factor;
36  int isp = species.nonadiabatic_idx;
37  Kokkos::parallel_for("add_f0_analytic", Kokkos::RangePolicy<HostExSpace>( 0, f0_f.n_nodes()), KOKKOS_LAMBDA(const int node){
38  // Add contribution from anomalous diffusion
39  int node_global = node+pol_decomp.node_offset;
40  double prefac2 = species.f0.temp_ev_h(node)/(species.f0.temp_ev_h(node) + species.f0.delta_T_h(node_global));
41  double prefac1 = sqrt(prefac2)*(species.f0.den_h(node) + species.f0.delta_n_h(node_global))
42  /(species.f0.temp_ev_h(node) + species.f0.delta_T_h(node_global));
43  double flow = species.f0.flow_h(node) + species.f0.delta_u_h(node_global)/species.get_f0_eq_thermal_velocity_lnode_h(node);
44 
45  double pot = 0.0;
46  if(species.is_electron){
47  // normalized potential
48 
49  // Get dpot
50  Simd<double> dpot;
51  const int i_simd = 0;
52  dpot.zero();
53 #ifdef XGCA
54 #ifdef COL_RELAX_TEST
55  // Keep dpot[i_simd] = 0
56 #else
57  // Upper limit for pot (We use linearized Poisson equation!!!)
58  // This is important at the wall strike points
59  electric_field.dpot_h(0, node_global).gather(dpot, i_simd, 1.0);
60 #endif
61 #elif defined(XGC1)
62  // Old Fortran remnant
64 
65  // Average the two planes
66  double wphi[2];
67  wphi[0] = 0.5;
68  wphi[1] = 0.5;
69  electric_field.dpot_ff_h(node_global).gather(dpot, i_simd, 1.0, wphi);
70 
71  // Old Fortran remnant
75 #endif
76  pot = dpot[i_simd]/(species.f0.temp_ev_h(node) + species.f0.delta_T_h(node_global));
77  pot = max(-sml.dpot_te_limit, min(sml.dpot_te_limit, pot));
78  }
79 
80  for (int ivr = 0; ivr < nvr; ivr++){
81  // normalized
82  double smu = ivr*dsmu; // mu or sqrt mu (normalized)
83  double mu = smu*smu; // now mu is normalized magnetic moment or (v_perp/v_th)^2
84 
85  for(int ivz = 0; ivz < nvz; ivz++){
86  // Shift ivz so it runs between -nvp and nvp
87  double vp = (ivz - nvp)*dvp - flow;
88 
89  if(ivr == 0){
90  smu = dsmu*inv_mu0_factor;
91  }
92  double en = 0.5*(vp*vp + mu);
93  // 'smu' is v_perp/v_th
94 
95  f0_f(isp, ivr, node, ivz) += sign*prefac1*smu*exp(-prefac2*en + pot);
96  }
97  }
98  });
99  Kokkos::fence();
101 }
102 
103 
114 
115  constexpr double eps=1e-2;
116 
117  int nvr = f0_f.n_vr(), nvz = f0_f.n_vz();
118  int nvp = (nvz - 1)/2;
119  double dsmu = f0_f.dsmu, dvp = f0_f.dvp;
120  double inv_mu0_factor = f0_f.inv_mu0_factor;
121 
123  int isp = species.nonadiabatic_idx;
124  Kokkos::parallel_for("f0_remove_negative", Kokkos::RangePolicy<HostExSpace>( 0, f0_f.n_nodes()), KOKKOS_LAMBDA(const int node){
125  double den_temp = species.f0.n_Ta_h(node);
126  double flow = species.f0.flow_h(node);
127  for (int ivr = 0; ivr < nvr; ivr++){
128  double smu = ivr*dsmu;
129  double mu = smu*smu;
130  if(ivr == 0){
131  smu = dsmu*inv_mu0_factor;
132  }
133  for(int ivz = 0; ivz < nvz; ivz++){
134  if(f0_f(isp, ivr, node, ivz) <= 0.0){
135  // Use <= to make sure we never have density = 0
136  // For the full distribution function, this
137  // is a small number
138  // Shift ivz so it runs between -nvp and nvp
139  double en = 0.5*( ((ivz - nvp)*dvp - flow)*
140  ((ivz - nvp)*dvp - flow) + mu);
141  f0_f(isp, ivr, node, ivz) = eps*smu*den_temp*exp(-en);
142  }
143  }
144  }
145  });
146  Kokkos::fence();
148 }
149 #endif
KOKKOS_INLINE_FUNCTION int n_vr() const
Definition: vgrid_distribution.hpp:153
Distribution< Device > f0
Species distribution in velocity space on local mesh nodes.
Definition: species.hpp:119
bool is_electron
Whether this species is the electrons.
Definition: species.hpp:78
Definition: sml.hpp:8
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)
Definition: add_f0_analytic.hpp:25
subroutine plasma(grid, itr, p, dene_out, deni_out, Te_out, Ti_out, Vparai_out)
Calculate the plasma density, temperature, and parallel velocity for a point in triangle itr using pl...
Definition: neutral_totalf.F90:1235
Definition: plasma.hpp:94
Definition: electric_field.hpp:36
void for_all_nonadiabatic_species(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:113
double dvp
grid spacing in parallel velocity
Definition: vgrid_distribution.hpp:25
int nonadiabatic_idx
Index of species skipping adiabatic species (for compatibility with fortran arrays) ...
Definition: species.hpp:80
int node_offset
Offset of first mesh node belonging to this MPI rank.
Definition: domain_decomposition.hpp:41
double inv_mu0_factor
Definition: vgrid_distribution.hpp:28
KOKKOS_INLINE_FUNCTION int n_vz() const
Definition: vgrid_distribution.hpp:161
KOKKOS_INLINE_FUNCTION void zero()
Definition: simd.hpp:76
double dpot_te_limit
Max absolute value of dpot/temp in getf0.
Definition: sml.hpp:96
double dsmu
grid spacing in mu
Definition: vgrid_distribution.hpp:27
KOKKOS_INLINE_FUNCTION double get_f0_eq_thermal_velocity_lnode_h(int inode) const
Definition: species.hpp:730
GridField< HostType, VarType::Scalar, PIT_GLOBAL, TorType::OnePlane, KinType::DriftKin > dpot_ff_h
Definition: electric_field.hpp:72
Definition: plasma.hpp:14
Definition: species.hpp:74
void parallel_for(const std::string name, int n_ptl, Function func, Option option, HostAoSoA aosoa_h, DeviceAoSoA aosoa_d)
Definition: streamed_parallel_for.hpp:252
KOKKOS_INLINE_FUNCTION int n_nodes() const
Definition: vgrid_distribution.hpp:157
GridField< HostType, VarType::Scalar, PhiInterpType::None, TorType::MultiplePlanes, KinType::DriftKin > dpot_h
Definition: electric_field.hpp:80
void f0_remove_negative(Plasma &plasma, VGridDistribution< HostType > &f0_f)
Definition: add_f0_analytic.hpp:113