10 extern "C" void init_ftn_lin(
int num,
double* psi,
double* var,
double psi_min,
double dpsi,
int n,
double* v);
28 template<
class Device>
30 static constexpr
int n = 1000;
34 Kokkos::View<double*,Kokkos::LayoutRight,Device>
v;
43 View<double*, CLayout, HostType> psi;
44 View<double*, CLayout, HostType> var;
48 file_reader.
open(filename);
49 file_reader.
read(&num, 1);
50 if(num<=0)
exit_XGC(
"\nError in profile ftn init : invalid number of data\n");
52 psi = View<double*, CLayout, HostType>(
NoInit(
"psi"), num);
53 var = View<double*, CLayout, HostType>(
NoInit(
"var"), num);
56 for(
int i=0; i<num; i++){
57 file_reader.
read(&psi(i), 1);
58 file_reader.
read(&var(i), 1);
63 file_reader.
read(&flag, 1);
64 if(flag!=-1)
exit_XGC(
"\nError in profile ftn init : invalid number of data or ending flag -1 is not set correctly\n");
73 psi = View<double*, CLayout, HostType>(
NoInit(
"psi"), num);
74 var = View<double*, CLayout, HostType>(
NoInit(
"var"), num);
78 MPI_Bcast( psi.data(), psi.size(), MPI_DOUBLE, root_rank,
SML_COMM_WORLD);
79 MPI_Bcast( var.data(), var.size(), MPI_DOUBLE, root_rank,
SML_COMM_WORLD);
83 for(
int i=0; i<num; i++){
84 psi(i) *= psi_normalization;
93 View<double*,CLayout,HostType> v_h(
NoInit(
"v_h"),
v.layout());
97 exit_XGC(
"\nError: CustomLinear profile shape requires PSPLINE.\n");
101 Kokkos::deep_copy(
v, v_h);
108 int num = var.size();
116 View<double*,CLayout,HostType> v_h(
NoInit(
"v_h"),
v.layout());
120 exit_XGC(
"\nError: CustomLinear profile shape requires PSPLINE.\n");
124 Kokkos::deep_copy(
v, v_h);
132 p_del((p_max_in-p_min_in)/(
n-1))
134 View<double*,CLayout,HostType> v_h(
NoInit(
"v_h"),
v.layout());
135 for (
int i = 0; i<
n; i++){
136 v_h(i) = 1.0 + i/(n-1.0);
140 Kokkos::deep_copy(
v, v_h);
145 KOKKOS_INLINE_FUNCTION
double interp(
double x)
const{
149 xn=max(min(xn,
double(
n)),0.0);
152 double wp=1.0-(xn-double(ip));
153 return v(ip)*wp+
v(ip+1)*(1.0-wp);
156 KOKKOS_INLINE_FUNCTION
double dinterp(
double x)
const{
159 xn=max(min(xn,
double(
n)),0.0);
170 template<
class Device>
173 static constexpr
int N_X = 5;
174 static constexpr
int N_Y = 4;
175 static constexpr
int N_V = 6;
192 return sv[1]*tanh((
inx[0]-psi)*
sv[2]) + sv[0];
194 return sv[0]*psi+
sv[1];
197 double tmp2=
sv[2]*(
inx[0]-sqrt(
inx[0]*psi));
198 double tmp3=exp(tmp2);
199 double tmp4=exp(-tmp2);
201 return sv[1]*((1+tmp2*
sv[3])*tmp3-tmp4)/(tmp3+tmp4) +
sv[0];
203 return lin.interp(psi);
216 }
else if(shape_int==1){
218 }
else if(shape_int==2){
220 }
else if(shape_int==5){
222 }
else if(shape_int==-1){
225 printf(
"Invalid shape number in profile setup: %d",shape_int);
226 exit_XGC(
"\nError: Invalid shape number\n");
231 Profile(
double psi_normalization,
Shape shape_in,
const std::vector<double>& inx_in,
const std::vector<double>& iny_in,
const std::string& filename)
235 for(
int i=0; i<inx_in.size(); i++)
inx[i] = inx_in[i];
236 for(
int i=0; i<iny_in.size(); i++)
iny[i] = iny_in[i];
260 /( - (
sv[1]+1.0e-99) *
sv[2] * sqrt(
inx[0]) );
264 exit_XGC(
"\nError: Invalid profile shape provided.\n");
269 Profile(
const View<double*, CLayout, HostType>& psi,
const View<double*, CLayout, HostType>& var)
278 for (
int i=0; i<5; i++)
inx[i]=0.0;
279 for (
int i=0; i<4; i++)
iny[i]=0.0;
280 for (
int i=0; i<6; i++)
sv[i]=0.0;
282 sv[1] = val_at_zero_psi;
289 for (
int i=0; i<5; i++)
inx[i]=0.0;
290 for (
int i=0; i<4; i++)
iny[i]=0.0;
291 for (
int i=0; i<6; i++)
sv[i]=0.0;
321 if(is_in_rgn12 && (psi_in<=b_field.
outpsi)){
323 }
else if(is_in_rgn12 && (psi_in>b_field.
outpsi)){
328 double psi=b_field.
outpsi;
332 }
else if(!is_in_rgn12 && z<b_field.
equil.
xpt.
z){
364 double tmp2 = 1.0/cosh((
inx[0]-psi_in)*
sv[2]);
365 return -sv[2]*sv[1]*tmp2*tmp2;
369 double tmp2=
sv[2]*(
inx[0]-sqrt(
inx[0]*psi_in));
370 double tmp3=exp(tmp2);
375 return -
inx[0]*
sv[1]*
sv[2]*tmp3*tmp3*(4.0+
sv[3]*(1.0+tmp3*tmp3+2.0*tmp2))
376 /(2.0*(1.0+tmp3*tmp3)*(1.0+tmp3*tmp3)*(sqrt(
inx[0]*psi_in)));
379 return lin.dinterp(psi_in);
385 void get_shape_inx_iny(
int& shape_int, std::vector<double>& inx_out, std::vector<double>& iny_out) {
389 for(
int i=0; i<
N_X; i++) inx_out[i] =
inx[i];
390 for(
int i=0; i<
N_Y; i++) iny_out[i] =
iny[i];
408 exit_XGC(
"\nError: Invalid profile shape provided.\n");
double p_del
Definition: profile.hpp:33
static constexpr int N_X
Definition: profile.hpp:173
bool is_rank_zero()
Definition: globals.hpp:27
MPI_Comm SML_COMM_WORLD
Definition: my_mpi.cpp:4
KOKKOS_INLINE_FUNCTION double slope(const MagneticField< DeviceType > &b_field, double psi_in, double r, double z) const
Definition: profile.hpp:352
KOKKOS_INLINE_FUNCTION double interp(double x) const
Definition: profile.hpp:145
double iny[N_Y]
Definition: profile.hpp:181
Shape shape
Definition: profile.hpp:179
RZPair xpt2
coordinates of 2nd X-point
Definition: equil.hpp:83
Definition: profile.hpp:17
static constexpr int N_V
Definition: profile.hpp:175
Definition: magnetic_field.hpp:14
double out_decay_width
width for exponential decay for psi>sml_outpsi
Definition: equil.hpp:91
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:46
Definition: profile.hpp:18
void read(T *var, int n)
Definition: file_reader.hpp:28
static Shape shape_from_int(int shape_int)
Definition: profile.hpp:213
CustomLinShape< Device > lin
Definition: profile.hpp:185
double z
Definition: grid_structs.hpp:30
bool set_xpt2
Whether to use a 2nd X-point.
Definition: equil.hpp:82
void get_shape_inx_iny(int &shape_int, std::vector< double > &inx_out, std::vector< double > &iny_out)
Definition: profile.hpp:385
Definition: profile.hpp:21
Definition: profile.hpp:20
double outpsi
Boundary condition used in a few spots.
Definition: magnetic_field.hpp:45
KOKKOS_INLINE_FUNCTION double value_at_coords(const MagneticField< DeviceType > &b_field, double psi, double r, double z) const
Definition: profile.hpp:188
double sv[N_V]
Definition: profile.hpp:182
double out_decay_factor
profiles decay exponentially to f(sml_outpsi)/decay_factor for psi>sml_outpsi
Definition: equil.hpp:89
Definition: profile.hpp:23
KOKKOS_INLINE_FUNCTION double dinterp(double x) const
Definition: profile.hpp:156
Definition: profile.hpp:19
Profile(double psi_normalization, Shape shape_in, const std::vector< double > &inx_in, const std::vector< double > &iny_in, const std::string &filename)
Definition: profile.hpp:231
Definition: profile.hpp:22
Profile()
Definition: profile.hpp:316
Definition: file_reader.hpp:6
double priv_flux_decay_factor
profiles decay exponentially to f(sml_outpsi)/decay_factor in priv. flux region
Definition: equil.hpp:90
double p_min
Definition: profile.hpp:31
int get_shape_int()
Definition: profile.hpp:396
static constexpr int N_Y
Definition: profile.hpp:174
double xpt2_psi
psi coordinate at 2nd X-point
Definition: equil.hpp:85
RZPair xpt
coordinates of 1st X-point
Definition: equil.hpp:80
CustomLinShape()
Definition: profile.hpp:143
double r
Definition: grid_structs.hpp:29
void exit_XGC(std::string msg)
Definition: globals.hpp:37
Shape
Definition: profile.hpp:15
void init_ftn_lin(int num, double *psi, double *var, double psi_min, double dpsi, int n, double *v)
Profile(double val_at_zero_psi, double slope)
Definition: profile.hpp:275
Definition: profile.hpp:24
KOKKOS_INLINE_FUNCTION double value(const MagneticField< DeviceType > &b_field, double psi_in, double r, double z) const
Definition: profile.hpp:319
Definition: profile.hpp:29
Profile(Shape shape_in)
Definition: profile.hpp:286
Profile(const View< double *, CLayout, HostType > &psi, const View< double *, CLayout, HostType > &var)
Definition: profile.hpp:269
double inx[N_X]
Definition: profile.hpp:180
KOKKOS_INLINE_FUNCTION bool is_in_region_1_or_2(double r, double z, double psi) const
Definition: magnetic_field.tpp:72
Definition: profile.hpp:16
static constexpr int n
Definition: profile.hpp:30
RZPair axis
Definition: equil.hpp:86
Kokkos::View< double *, Kokkos::LayoutRight, Device > v
Definition: profile.hpp:34
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:68
Definition: profile.hpp:171
double p_max
Definition: profile.hpp:32
double priv_flux_decay_width
width for exponential decay in private flux region
Definition: equil.hpp:92
void open(const std::string &input)
Definition: file_reader.hpp:13
double xpt_psi
Psi coordinate of 1st X-point.
Definition: equil.hpp:78