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 or no poloidal decomposition
40  den_all(NoInit("den_all"), (do_gather && pol_decomp.pol_decomp) ? grid.nnode : 1, 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){
51  if(pol_decomp.pol_decomp){
52  gather(pol_decomp);
53  }else{
54  // No gather is needed if there is no domain decomposition. Instead, just point the _all views to the local views
60  }
61  }
62 
63  // Set fortran pointers
64  set_mom_module_ptrs(den_local.data(), upara_local.data(), uperp_local.data(), temp_local.data(), Kperp_local.data(),
65  Kpara_local.data(), temp_para_local.data(), den_all.data(), temp_all.data(), temp_perp_all.data(),
66  temp_para_all.data(), upara_all.data());
67  }
68 
70  // Create local copies of these member views so that they can be passed into the lambda function
71  auto den = den_local;
72  auto upara = upara_local;
73  auto uperp = uperp_local;
74  auto temp = temp_local;
75  auto Kperp = Kperp_local;
76  auto Kpara = Kpara_local;
77  auto temp_para = temp_para_local;
78 
80  int isp = species.nonadiabatic_idx;
81 
82  // Calculate sum of f, E*f, vp*f
83  Kokkos::parallel_for("calculate_f0_charge_den", Kokkos::RangePolicy<HostExSpace>( 0, f.n_nodes()), KOKKOS_LAMBDA(const int inode){
84  for (int ivr=0; ivr<vgrid.nvr; ivr++){
85  // normalized volume element
86  double mu_vol = vgrid.mu_vol_fac(ivr);
87  double smu = vgrid.vperp_n(ivr);
88  double mu = smu*smu; // <--- This is (v_perp/v_th)^2 for v_perp grid and mu_N for sqrt(mu) grid!!!!
89  double en_perp= 0.5 * mu; // energy normalized by T, 0.5*(vperp/vth)^2
90 
91  // Sum over parallel velocity
92  for(int ivz=0; ivz<vgrid.nvz; ivz++){
93  double vp = vgrid.vpar_n(ivz);
94  double en_para2 = vp*vp; // energy normalized by T, 0.5*(vpara/vth)^2
95 
96  double wt = mu_vol*f(isp, ivr, inode, ivz);
97  den(inode,isp) += wt;
98  upara(inode,isp) += wt*vp;
99  uperp(inode,isp) += wt*smu;
100  Kperp(inode,isp) += wt*en_perp;
101  Kpara(inode,isp) += wt*en_para2;
102  }
103  }
104  });
105  Kokkos::fence();
106 
107  // Normalization of all moments and compute temp_para and temp
108  const double UPERP_DEFAULT_COEFF = sqrt(0.25*TWOPI);
109  Kokkos::parallel_for("calculate_f0_charge_den", Kokkos::RangePolicy<HostExSpace>( 0, f.n_nodes()), KOKKOS_LAMBDA(const int inode){
110 
111  // Impose lower limit on density
112  double den0 = 1.0e-4*species.f0.den_h(inode); // For safety, to prevent inf and NaN
113  if (!std::isfinite(den(inode,isp))) den(inode,isp) = den0;
114  den(inode,isp) = max(den0,den(inode,isp));
115 
116  // Normalized v
117  double t_norm_ev = species.f0.fg_temp_ev_h(inode);
118  double vth = 1./species.f0.fg_vth_inv_h(inode);
119 
120  // Normalize by density
121  double inv_density = 1.0/den(inode,isp);
122  Kperp(inode,isp) *= t_norm_ev*inv_density;
123  Kpara(inode,isp) *= t_norm_ev*inv_density;
124  upara(inode,isp) *= vth*inv_density;
125  uperp(inode,isp) *= vth*inv_density;
126  den(inode,isp) *= species.f0.grid_vol_vonly_h(inode);
127 
128  // Calculate temp_para and temp
129  temp_para(inode,isp) = Kpara(inode,isp) - 2.0*kinetic_energy(species.mass, upara(inode,isp))*J_2_EV;
130  temp(inode,isp) = ( temp_para(inode,isp) + 2.0*Kperp(inode,isp) )/3.0;
131 
132  // Set possible NaNs and Infs to default values
133  if (!std::isfinite(den(inode,isp))) den(inode,isp) = species.f0.den_h(inode);
134  if (!std::isfinite(Kperp(inode,isp))) Kperp(inode,isp) = t_norm_ev;
135  if (!std::isfinite(Kpara(inode,isp))) Kpara(inode,isp) = t_norm_ev;
136  if (!std::isfinite(temp_para(inode,isp))) temp_para(inode,isp) = t_norm_ev;
137  if (!std::isfinite(temp(inode,isp))) temp(inode,isp) = t_norm_ev;
138  if (!std::isfinite(upara(inode,isp))) upara(inode,isp) = species.f0.flow_h(inode);
139  if (!std::isfinite(uperp(inode,isp))) uperp(inode,isp) = UPERP_DEFAULT_COEFF*vth;
140 
141  // Impose lower and upper limits
142  den (inode,isp) = min(max(den (inode,isp),1.0e-2*species.f0.den_h(inode)),1.0e2*species.f0.den_h(inode));
143  Kperp (inode,isp) = min(max(Kperp (inode,isp),1.0e-2*t_norm_ev ),1.0e2*t_norm_ev );
144  Kpara (inode,isp) = min(max(Kpara (inode,isp),1.0e-2*t_norm_ev ),1.0e2*t_norm_ev );
145  temp (inode,isp) = min(max(temp (inode,isp),1.0e-2*t_norm_ev ),1.0e2*t_norm_ev );
146  temp_para (inode,isp) = min(max(temp_para (inode,isp),1.0e-2*t_norm_ev ),1.0e2*t_norm_ev );
147  upara (inode,isp) = min(max(upara (inode,isp),-0.75*vgrid.vp_max*vth ),0.75*vgrid.vp_max*vth );
148  uperp (inode,isp) = min(max(uperp (inode,isp),0.05*vgrid.vp_max*vth ),0.75*vgrid.vp_max*vth );
149  });
150  Kokkos::fence();
152  }
153 
154  // Gathers all moments on local nodes, and distributes copies to all compute nodes
155  void gather(const DomainDecomposition<DeviceType>& pol_decomp){
156 #ifdef USE_MPI
157  View<int*,CLayout,HostType> cnts(NoInit("cnts"), pol_decomp.mpi.n_plane_ranks);
158  View<int*,CLayout,HostType> displs(NoInit("cnts"), pol_decomp.mpi.n_plane_ranks);
159 
160  // Create cnts array with # vertex nodes on each processor
161  int nsp = den_local.extent(1);
162  for(int i = 0; i<cnts.size(); i++){
163  displs(i) = nsp*(pol_decomp.gvid0_pid_h(i) - 1); // gvid0_pid_h is 1-indexed
164  cnts(i) = nsp*(pol_decomp.gvid0_pid_h(i+1) - pol_decomp.gvid0_pid_h(i));
165  }
166 
167  int my_send_count = cnts(pol_decomp.mpi.my_plane_rank);
168  MPI_Allgatherv(den_local.data(), my_send_count, MPI_DOUBLE, den_all.data(), cnts.data(), displs.data(), MPI_DOUBLE, pol_decomp.mpi.plane_comm);
169  MPI_Allgatherv(temp_local.data(), my_send_count, MPI_DOUBLE, temp_all.data(), cnts.data(), displs.data(), MPI_DOUBLE, pol_decomp.mpi.plane_comm);
170  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);
171  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);
172  MPI_Allgatherv(upara_local.data(), my_send_count, MPI_DOUBLE, upara_all.data(), cnts.data(), displs.data(), MPI_DOUBLE, pol_decomp.mpi.plane_comm);
173 #endif
174  }
175 };
176 
177 #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:69
Moments()
Definition: moments.hpp:28
Distribution< Device > f0
Species distribution in velocity space on local mesh nodes.
Definition: species.hpp:127
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:1224
Definition: plasma.hpp:113
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:132
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 mu_vol_fac(int ivr) const
Definition: velocity_grid.hpp:80
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
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:155
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:159
View< double **, CLayout, HostType > uperp_local
perpendicular bulk velocity, local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:16
bool pol_decomp
Use poloidal decomposition.
Definition: domain_decomposition.hpp:46
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:63