XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
moments.hpp
Go to the documentation of this file.
1 #ifndef MOMENTS_HPP
2 #define MOMENTS_HPP
3 
4 #include "grid.hpp"
6 #include "vgrid_distribution.hpp"
7 #include "plasma.hpp"
8 #include "basic_physics.hpp"
9 
10 extern "C" void set_mom_module_ptrs(double* den_local_cpp, double* upara_local_cpp, double* uperp_local_cpp, double* temp_local_cpp, double* Kperp_local_cpp,
11  double* Kpara_local_cpp, double* temp_para_local_cpp, double* den_all_cpp, double* temp_all_cpp, double* temp_perp_all_cpp, double* temp_para_all_cpp, double* upara_all_cpp);
12 
13 struct Moments{
14  View<double**, CLayout, HostType> den_local;
15  View<double**, CLayout, HostType> upara_local;
16  View<double**, CLayout, HostType> uperp_local;
17  View<double**, CLayout, HostType> temp_local;
18  View<double**, CLayout, HostType> Kperp_local;
19  View<double**, CLayout, HostType> Kpara_local;
20  View<double**, CLayout, HostType> temp_para_local;
21 
22  View<double**, CLayout, HostType> den_all;
23  View<double**, CLayout, HostType> upara_all;
24  View<double**, CLayout, HostType> temp_all;
25  View<double**, CLayout, HostType> temp_perp_all;
26  View<double**, CLayout, HostType> temp_para_all;
27 
28  Moments(){}
29 
30  Moments(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const VelocityGrid& vgrid, Plasma& plasma, const VGridDistribution<HostType>& f, bool do_gather)
31  : den_local("den_local", f.n_nodes(), f.n_species()),
32  upara_local("upara_local", den_local.layout()),
33  uperp_local("uperp_local", den_local.layout()),
34  temp_local(NoInit("temp_local"), den_local.layout()), // NoInit because it is derived from the other quantities
35  Kperp_local("Kperp_local", den_local.layout()),
36  Kpara_local("Kpara_local", den_local.layout()),
37  temp_para_local(NoInit("temp_para_local"), den_local.layout()), // NoInit because it is derived from the other quantities
38 
39  // These don't need to be allocated if do_gather==false, but allocate anyway for now for Fortran pointers
40  den_all(NoInit("den_all"), grid.nnode, f.n_species()),
41  upara_all(NoInit("upara_all"), den_all.layout()),
42  temp_all(NoInit("temp_all"), den_all.layout()),
43  temp_perp_all(NoInit("temp_perp_all"), den_all.layout()),
44  temp_para_all(NoInit("temp_para_all"), den_all.layout())
45  {
46  // Calculate moments of f in local domain
47  calculate_moments_of_distribution(vgrid, plasma, f);
48 
49  // Gather decomposed moments if requested
50  if(do_gather) gather(pol_decomp);
51 
52  // Set fortran pointers
53  set_mom_module_ptrs(den_local.data(), upara_local.data(), uperp_local.data(), temp_local.data(), Kperp_local.data(),
54  Kpara_local.data(), temp_para_local.data(), den_all.data(), temp_all.data(), temp_perp_all.data(),
55  temp_para_all.data(), upara_all.data());
56  }
57 
59  // Create local copies of these member views so that they can be passed into the lambda function
60  auto den = den_local;
61  auto upara = upara_local;
62  auto uperp = uperp_local;
63  auto temp = temp_local;
64  auto Kperp = Kperp_local;
65  auto Kpara = Kpara_local;
66  auto temp_para = temp_para_local;
67 
69  int isp = species.nonadiabatic_idx;
70 
71  // Calculate sum of f, E*f, vp*f
72  Kokkos::parallel_for("calculate_f0_charge_den", Kokkos::RangePolicy<HostExSpace>( 0, f.n_nodes()), KOKKOS_LAMBDA(const int inode){
73  for (int ivr=0; ivr<vgrid.nvr; ivr++){
74  // normalized volume element
75  double mu_vol = (ivr==0 || ivr==vgrid.nvr-1) ? 0.5 : 1.0;
76  double smu = vgrid.vperp_n(ivr);
77  double mu = smu*smu; // <--- This is (v_perp/v_th)^2 for v_perp grid and mu_N for sqrt(mu) grid!!!!
78  double en_perp= 0.5 * mu; // energy normalized by T, 0.5*(vperp/vth)^2
79 
80  // Sum over parallel velocity
81  for(int ivz=0; ivz<vgrid.nvz; ivz++){
82  double vp = vgrid.vpar_n(ivz);
83  double en_para2 = vp*vp; // energy normalized by T, 0.5*(vpara/vth)^2
84 
85  double wt = mu_vol*f(isp, ivr, inode, ivz);
86  den(inode,isp) += wt;
87  upara(inode,isp) += wt*vp;
88  uperp(inode,isp) += wt*smu;
89  Kperp(inode,isp) += wt*en_perp;
90  Kpara(inode,isp) += wt*en_para2;
91  }
92  }
93  });
94  Kokkos::fence();
95 
96  // Normalization of all moments and compute temp_para and temp
97  const double UPERP_DEFAULT_COEFF = sqrt(0.25*TWOPI);
98  Kokkos::parallel_for("calculate_f0_charge_den", Kokkos::RangePolicy<HostExSpace>( 0, f.n_nodes()), KOKKOS_LAMBDA(const int inode){
99 
100  // Impose lower limit on density
101  double den0 = 1.0e-4*species.f0.den_h(inode); // For safety, to prevent inf and NaN
102  if (!std::isfinite(den(inode,isp))) den(inode,isp) = den0;
103  den(inode,isp) = max(den0,den(inode,isp));
104 
105  // Normalized v
106  double t_norm_ev = species.f0.temp_ev_h(inode);
107  double vth = species.get_f0_eq_thermal_velocity_lnode_h(inode);
108 
109  // Normalize by density
110  double inv_density = 1.0/den(inode,isp);
111  Kperp(inode,isp) *= t_norm_ev*inv_density;
112  Kpara(inode,isp) *= t_norm_ev*inv_density;
113  upara(inode,isp) *= vth*inv_density;
114  uperp(inode,isp) *= vth*inv_density;
115  den(inode,isp) *= species.f0.grid_vol_vonly_h(inode);
116 
117  // Calculate temp_para and temp
118  temp_para(inode,isp) = Kpara(inode,isp) - 2.0*kinetic_energy(species.mass, upara(inode,isp))*J_2_EV;
119  temp(inode,isp) = ( temp_para(inode,isp) + 2.0*Kperp(inode,isp) )/3.0;
120 
121  // Set possible NaNs and Infs to default values
122  if (!std::isfinite(den(inode,isp))) den(inode,isp) = species.f0.den_h(inode);
123  if (!std::isfinite(Kperp(inode,isp))) Kperp(inode,isp) = t_norm_ev;
124  if (!std::isfinite(Kpara(inode,isp))) Kpara(inode,isp) = t_norm_ev;
125  if (!std::isfinite(temp_para(inode,isp))) temp_para(inode,isp) = t_norm_ev;
126  if (!std::isfinite(temp(inode,isp))) temp(inode,isp) = t_norm_ev;
127  if (!std::isfinite(upara(inode,isp))) upara(inode,isp) = species.f0.flow_h(inode)*vth;
128  if (!std::isfinite(uperp(inode,isp))) uperp(inode,isp) = UPERP_DEFAULT_COEFF*vth;
129 
130  // Impose lower and upper limits
131  den (inode,isp) = min(max(den (inode,isp),1.0e-2*species.f0.den_h(inode)),1.0e2*species.f0.den_h(inode));
132  Kperp (inode,isp) = min(max(Kperp (inode,isp),1.0e-2*t_norm_ev ),1.0e2*t_norm_ev );
133  Kpara (inode,isp) = min(max(Kpara (inode,isp),1.0e-2*t_norm_ev ),1.0e2*t_norm_ev );
134  temp (inode,isp) = min(max(temp (inode,isp),1.0e-2*t_norm_ev ),1.0e2*t_norm_ev );
135  temp_para (inode,isp) = min(max(temp_para (inode,isp),1.0e-2*t_norm_ev ),1.0e2*t_norm_ev );
136  upara (inode,isp) = min(max(upara (inode,isp),-0.75*vgrid.vp_max*vth ),0.75*vgrid.vp_max*vth );
137  uperp (inode,isp) = min(max(uperp (inode,isp),0.05*vgrid.vp_max*vth ),0.75*vgrid.vp_max*vth );
138  });
139  Kokkos::fence();
141  }
142 
143  // Gathers all moments on local nodes, and distributes copies to all compute nodes
144  void gather(const DomainDecomposition<DeviceType>& pol_decomp){
145 #ifdef USE_MPI
146  View<int*,CLayout,HostType> cnts(NoInit("cnts"), pol_decomp.mpi.n_plane_ranks);
147  View<int*,CLayout,HostType> displs(NoInit("cnts"), pol_decomp.mpi.n_plane_ranks);
148 
149  // Create cnts array with # vertex nodes on each processor
150  int nsp = den_local.extent(1);
151  for(int i = 0; i<cnts.size(); i++){
152  displs(i) = nsp*(pol_decomp.gvid0_pid_h(i) - 1); // gvid0_pid_h is 1-indexed
153  cnts(i) = nsp*(pol_decomp.gvid0_pid_h(i+1) - pol_decomp.gvid0_pid_h(i));
154  }
155 
156  int my_send_count = cnts(pol_decomp.mpi.my_plane_rank);
157  MPI_Allgatherv(den_local.data(), my_send_count, MPI_DOUBLE, den_all.data(), cnts.data(), displs.data(), MPI_DOUBLE, pol_decomp.mpi.plane_comm);
158  MPI_Allgatherv(temp_local.data(), my_send_count, MPI_DOUBLE, temp_all.data(), cnts.data(), displs.data(), MPI_DOUBLE, pol_decomp.mpi.plane_comm);
159  MPI_Allgatherv(temp_para_local.data(),my_send_count, MPI_DOUBLE, temp_para_all.data(), cnts.data(), displs.data(), MPI_DOUBLE, pol_decomp.mpi.plane_comm);
160  MPI_Allgatherv(Kperp_local.data(), my_send_count, MPI_DOUBLE, temp_perp_all.data(), cnts.data(), displs.data(), MPI_DOUBLE, pol_decomp.mpi.plane_comm);
161  MPI_Allgatherv(upara_local.data(), my_send_count, MPI_DOUBLE, upara_all.data(), cnts.data(), displs.data(), MPI_DOUBLE, pol_decomp.mpi.plane_comm);
162 #endif
163  }
164 };
165 
166 #endif
View< double **, CLayout, HostType > Kperp_local
perpendicular kinetic energy, , local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:18
void calculate_moments_of_distribution(const VelocityGrid &vgrid, Plasma &plasma, const VGridDistribution< HostType > &f)
Definition: moments.hpp:58
Moments()
Definition: moments.hpp:28
Distribution< Device > f0
Species distribution in velocity space on local mesh nodes.
Definition: species.hpp:120
View< double **, CLayout, HostType > temp_local
temperature, , local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:17
View< double **, CLayout, HostType > temp_all
temperature, , on all nodes [Nsp,Nnode]
Definition: moments.hpp:24
Definition: velocity_grid.hpp:8
View< double **, CLayout, HostType > upara_all
parallel bulk velocity, , on all nodes [Nsp,Nnode]
Definition: moments.hpp:23
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:1237
Definition: plasma.hpp:109
Definition: moments.hpp:13
int nvr
full grid size (including zero)
Definition: velocity_grid.hpp:19
void for_all_nonadiabatic_species(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:128
View< double **, CLayout, HostType > den_all
density, , on all nodes [Nsp,Nnode]
Definition: moments.hpp:22
int nonadiabatic_idx
Index of species skipping adiabatic species (for compatibility with fortran arrays) ...
Definition: species.hpp:81
KOKKOS_INLINE_FUNCTION double vpar_n(int ivz) const
Definition: velocity_grid.hpp:70
View< double **, CLayout, HostType > upara_local
parallel bulk velocity, local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:15
View< double **, CLayout, HostType > Kpara_local
parllel kinetic energy, , local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:19
constexpr double J_2_EV
Conversion rate J to ev.
Definition: constants.hpp:6
double mass
Particle mass.
Definition: species.hpp:83
KOKKOS_INLINE_FUNCTION double kinetic_energy(double mass, double v)
Definition: basic_physics.hpp:8
View< double **, CLayout, HostType > temp_perp_all
perpendicular kinetic energy, , on all nodes [Nsp,Nnode]
Definition: moments.hpp:25
View< double **, CLayout, HostType > den_local
density, local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:14
double vp_max
max parallel velocity
Definition: velocity_grid.hpp:10
KOKKOS_INLINE_FUNCTION double get_f0_eq_thermal_velocity_lnode_h(int inode) const
Definition: species.hpp:727
Definition: plasma.hpp:14
View< double **, CLayout, HostType > temp_para_local
parallel temperature, , local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:20
void gather(const DomainDecomposition< DeviceType > &pol_decomp)
Definition: moments.hpp:144
void set_mom_module_ptrs(double *den_local_cpp, double *upara_local_cpp, double *uperp_local_cpp, double *temp_local_cpp, double *Kperp_local_cpp, double *Kpara_local_cpp, double *temp_para_local_cpp, double *den_all_cpp, double *temp_all_cpp, double *temp_perp_all_cpp, double *temp_para_all_cpp, double *upara_all_cpp)
View< double **, CLayout, HostType > temp_para_all
parllel kinetic energy, , on all nodes [Nsp,Nnode]
Definition: moments.hpp:26
Definition: species.hpp:75
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
int nvz
full grid size (including negative and zero)
Definition: velocity_grid.hpp:20
KOKKOS_INLINE_FUNCTION int n_nodes() const
Definition: vgrid_distribution.hpp:157
View< double **, CLayout, HostType > uperp_local
perpendicular bulk velocity, local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:16
Moments(const Grid< DeviceType > &grid, const DomainDecomposition< DeviceType > &pol_decomp, const VelocityGrid &vgrid, Plasma &plasma, const VGridDistribution< HostType > &f, bool do_gather)
Definition: moments.hpp:30
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:68
constexpr double TWOPI
Definition: constants.hpp:9
KOKKOS_INLINE_FUNCTION double vperp_n(int ivr) const
Definition: velocity_grid.hpp:75
Kokkos::View< int *, Kokkos::LayoutRight, HostType > gvid0_pid_h
Which processors get which vertices (host)
Definition: domain_decomposition.hpp:59