1 #ifndef ANALYTIC_PSI_HPP
2 #define ANALYTIC_PSI_HPP
12 double get_circular_analytic_psi(
const Equilibrium& equil,
double r,
double z){
20 double get_analytic_psi(
const Equilibrium& equil,
double r,
double z){
21 double rp = r - equil.
axis.
r;
22 double zp = z - equil.
axis.
z;
23 double r_xpt1p = equil.
xpt.
r-equil.
axis.
r;
24 double z_xpt1p = equil.
xpt.
z-equil.
axis.
z;
27 double zx1p = z - equil.
xpt.
z;
28 double zx2p = z - equil.
xpt2.
z;
29 double avg_z_xpt = (equil.
xpt2.
z - equil.
xpt.
z)/2;
31 double ar = 1/pow(avg_z_xpt,2);
33 double r_xptp = (zp>0 ? r_xpt2p : r_xpt1p);
34 double z_xptp = (zp>0 ? z_xpt2p : z_xpt1p);
35 double zxp = (zp>0 ? zx2p : zx1p);
38 double a = -2/pow(z_xptp,3);
39 double b = -3/pow(z_xptp,2);
40 if (zp>z_xpt2p || zp<z_xpt1p){
47 double no_xpt_psi = 0;
53 anx = 1/pow(z_xpt1p,2);
59 return ar*pow(rp-r_xptp/z_xptp*zp,2)*avg_psi_xpt + (a*pow(zxp,3) + b*pow(zxp,2))*xpt_psi + xpt_psi
60 + anx*pow(zp,2)*no_xpt_psi - no_xpt_psi;
63 void get_analytic_psi_derivs(
const Equilibrium& equil,
double r,
double z,
double& dpsidr,
double& dpsidz){
64 double rp = r - equil.
axis.
r;
65 double zp = z - equil.
axis.
z;
66 double r_xpt1p = equil.
xpt.
r-equil.
axis.
r;
67 double z_xpt1p = equil.
xpt.
z-equil.
axis.
z;
70 double zx1p = z - equil.
xpt.
z;
71 double zx2p = z - equil.
xpt2.
z;
72 double avg_z_xpt = (equil.
xpt2.
z - equil.
xpt.
z)/2;
74 double ar = 1/pow(avg_z_xpt,2);
76 double r_xptp = (zp>0 ? r_xpt2p : r_xpt1p);
77 double z_xptp = (zp>0 ? z_xpt2p : z_xpt1p);
78 double zxp = (zp>0 ? zx2p : zx1p);
81 double a = -2/pow(z_xptp,3);
82 double b = -3/pow(z_xptp,2);
83 if (zp>z_xpt2p || zp<z_xpt1p){
90 double no_xpt_psi = 0;
96 anx = 1/pow(z_xpt1p,2);
102 dpsidr = 2*ar*(rp-r_xptp/z_xptp*zp)*avg_psi_xpt;
103 dpsidz = -2*ar*(rp-r_xptp/z_xptp*zp)*avg_psi_xpt*r_xptp/z_xptp + (3*a*pow(zxp,2) + 2*b*zxp)*xpt_psi + 2*anx*zp*no_xpt_psi;
107 void get_analytic_psi_second_derivs(
const Equilibrium& equil,
double z,
double& dpsidrdz,
double& dpsidr2,
double& dpsidz2){
108 double zp = z - equil.
axis.
z;
109 double r_xpt1p = equil.
xpt.
r-equil.
axis.
r;
110 double z_xpt1p = equil.
xpt.
z-equil.
axis.
z;
113 double zx1p = z - equil.
xpt.
z;
114 double zx2p = z - equil.
xpt2.
z;
115 double avg_z_xpt = (equil.
xpt2.
z - equil.
xpt.
z)/2;
117 double ar = 1/pow(avg_z_xpt,2);
119 double r_xptp = (zp>0 ? r_xpt2p : r_xpt1p);
120 double z_xptp = (zp>0 ? z_xpt2p : z_xpt1p);
121 double zxp = (zp>0 ? zx2p : zx1p);
124 double a = -2/pow(z_xptp,3);
125 double b = -3/pow(z_xptp,2);
126 if (zp>z_xpt2p || zp<z_xpt1p){
128 b = -1/pow(z_xptp,2);
133 double no_xpt_psi = 0;
139 anx = 1/pow(z_xpt1p,2);
145 dpsidrdz = -2*ar*(r_xptp/z_xptp)*avg_psi_xpt;
146 dpsidr2 = 2*ar*avg_psi_xpt;
147 dpsidz2 = 2*ar*(r_xptp/z_xptp)*avg_psi_xpt*r_xptp/z_xptp + (6*a*zxp + 2*b)*xpt_psi + 2*anx*no_xpt_psi;
156 std::vector<int> prevrow(nx+1,-1);
158 for (
int j=ny;j>=0;j--){
161 for (
int i=0;i<=nx;i++){
164 int num = get_analytic_psi(equil,r,z)*5.0;
165 char mychar = 48+num;
166 if (prevnum!=num || prevrow[i]!=num) printf(
"%c",mychar);
176 Kokkos::View<double**, Kokkos::LayoutRight, HostType> analytic_psi_grid(
const Equilibrium& equil,
177 const View<double*, HostType>& r,
178 const View<double*, HostType>& z){
180 Kokkos::View<double**, Kokkos::LayoutRight, HostType> psi_grid(
"psi_grid",z.size(),r.size());
181 for (
int j=0; j<z.size(); j++){
182 for (
int i=0; i<r.size(); i++){
183 psi_grid(j,i) = get_analytic_psi(equil, r(i), z(j));
190 Kokkos::View<double**, Kokkos::LayoutRight, HostType> circular_analytic_psi_grid(
const Equilibrium& equil,
191 const View<double*, HostType>& r,
192 const View<double*, HostType>& z){
194 Kokkos::View<double**, Kokkos::LayoutRight, HostType> psi_grid(
"psi_grid",z.size(),r.size());
195 for (
int j=0; j<z.size(); j++){
196 for (
int i=0; i<r.size(); i++){
197 psi_grid(j,i) = get_circular_analytic_psi(equil, r(i), z(j));
RZPair xpt2
coordinates of 2nd X-point
Definition: equil.hpp:83
double z
Definition: grid_structs.hpp:30
bool set_xpt2
Whether to use a 2nd X-point.
Definition: equil.hpp:82
double min_r
Definition: rz_bounds.hpp:5
RZBounds bounds
Min and max for r and z.
Definition: equil.hpp:77
double xpt2_psi
psi coordinate at 2nd X-point
Definition: equil.hpp:85
RZPair xpt
coordinates of 1st X-point
Definition: equil.hpp:80
double r
Definition: grid_structs.hpp:29
KOKKOS_INLINE_FUNCTION double get_z_width() const
Definition: rz_bounds.hpp:28
KOKKOS_INLINE_FUNCTION double get_r_width() const
Definition: rz_bounds.hpp:24
double min_z
Definition: rz_bounds.hpp:7
RZPair axis
Definition: equil.hpp:86
double xpt_psi
Psi coordinate of 1st X-point.
Definition: equil.hpp:78