Skip to content

Commit b738f5a

Browse files
cianciosacianciosa
cianciosa
authored and
cianciosa
committed
Add explicit O and X mode dipersion relations.
1 parent 54c2d0d commit b738f5a

File tree

3 files changed

+144
-2
lines changed

3 files changed

+144
-2
lines changed

‎graph_driver/xrays.cpp‎

Lines changed: 25 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -506,6 +506,28 @@ void trace_ray(const commandline::parser &cl,
506506
stream.str(),
507507
local_num_rays,
508508
thread_number);
509+
} else if (dispersion == "cold_plasma-O") {
510+
run_dispersion<dispersion::single_mode_wave<T, SAFE_MATH, 1>> (cl, omega,
511+
kx, ky, kz,
512+
x, y, z,
513+
t, dt, eq,
514+
num_steps,
515+
sub_steps,
516+
engine,
517+
stream.str(),
518+
local_num_rays,
519+
thread_number);
520+
} else if (dispersion == "cold_plasma-X") {
521+
run_dispersion<dispersion::single_mode_wave<T, SAFE_MATH, -1>> (cl, omega,
522+
kx, ky, kz,
523+
x, y, z,
524+
t, dt, eq,
525+
num_steps,
526+
sub_steps,
527+
engine,
528+
stream.str(),
529+
local_num_rays,
530+
thread_number);
509531
} else {
510532
run_dispersion<dispersion::cold_plasma<T, SAFE_MATH>> (cl, omega,
511533
kx, ky, kz,
@@ -972,7 +994,9 @@ commandline::parser parse_commandline(int argc, const char * argv[]) {
972994
"bohm_gross",
973995
"ordinary_wave",
974996
"extra_ordinary_wave",
975-
"cold_plasma"
997+
"cold_plasma",
998+
"cold_plasma-O",
999+
"cold_plasma-X"
9761000
});
9771001
cl.add_option("equilibrium", true, "Equilibrium to use.", {
9781002
"efit",

‎graph_framework/dispersion.hpp‎

Lines changed: 118 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -877,6 +877,124 @@ namespace dispersion {
877877
}
878878
};
879879

880+
//------------------------------------------------------------------------------
881+
/// @brief Ordinary wave dispersion function.
882+
///
883+
/// @tparam T Base type of the calculation.
884+
/// @tparam SAFE_MATH Use safe math operations.
885+
/// @tparam MODE Mode selector -1 X-mode, +1 O-mode.
886+
//------------------------------------------------------------------------------
887+
template<jit::float_scalar T, bool SAFE_MATH=false, int MODE=1>
888+
class single_mode_wave final : public physics<T, SAFE_MATH> {
889+
public:
890+
//------------------------------------------------------------------------------
891+
/// @brief Cold plasma dispersion relation for a single mode.
892+
///
893+
/// D = N^2 - (-B + MODE*Sqrt(B^2 - 4AC))/(2A) (1)
894+
///
895+
/// A = ε_1*sin^2(θ) + ε_3*cos^2(θ) (2)
896+
///
897+
/// B = -ε_1*ε_3*(1 + cos^2(θ)) - (ε_1^2 - ε_2^2)*sin^2(θ) (3)
898+
///
899+
/// C = ε_3*(ε_1^2 + ε_2^2) (4)
900+
///
901+
/// θ is defined as
902+
/// cos(θ) = B·N/BN (5)
903+
///
904+
/// ε_1 = ε⟂ = 1 - Sum_i=1 X_i/(1 - Y^2_i) (6)
905+
///
906+
/// ε_2 = g = -X_eY_e/(1 - Y_e^2) + Sum_i=2 X_iY_i/(1 - Y^2_i) (7)
907+
///
908+
/// ε_3 = ε∥ = 1 - Sum_i=1 X_i (8)
909+
///
910+
/// X_i = ωpi^2/ω^2 (9)
911+
///
912+
/// Y_i = ωci/ω (10)
913+
///
914+
/// ωpi is the plasma frequency and ωci is the cyclotron frequency.
915+
///
916+
/// @params[in] w Omega variable.
917+
/// @params[in] kx Kx variable.
918+
/// @params[in] ky Ky variable.
919+
/// @params[in] kz Kz variable.
920+
/// @params[in] x x variable.
921+
/// @params[in] y y variable.
922+
/// @params[in] z z variable.
923+
/// @params[in] t Current time.
924+
/// @params[in] eq The plasma equilibrium.
925+
//------------------------------------------------------------------------------
926+
virtual graph::shared_leaf<T, SAFE_MATH>
927+
D(graph::shared_leaf<T, SAFE_MATH> w,
928+
graph::shared_leaf<T, SAFE_MATH> kx,
929+
graph::shared_leaf<T, SAFE_MATH> ky,
930+
graph::shared_leaf<T, SAFE_MATH> kz,
931+
graph::shared_leaf<T, SAFE_MATH> x,
932+
graph::shared_leaf<T, SAFE_MATH> y,
933+
graph::shared_leaf<T, SAFE_MATH> z,
934+
graph::shared_leaf<T, SAFE_MATH> t,
935+
equilibrium::shared<T, SAFE_MATH> &eq) {
936+
// Dielectric terms.
937+
// Frequencies
938+
auto ne = eq->get_electron_density(x, y, z);
939+
auto wpe2 = build_plasma_fequency(ne,
940+
physics<T, SAFE_MATH>::q,
941+
physics<T, SAFE_MATH>::me,
942+
physics<T, SAFE_MATH>::c,
943+
physics<T, SAFE_MATH>::epsion0);
944+
auto b_vec = eq->get_magnetic_field(x, y, z);
945+
auto b_len = b_vec->length();
946+
auto ec = build_cyclotron_fequency(-physics<T, SAFE_MATH>::q,
947+
b_len,
948+
physics<T, SAFE_MATH>::me,
949+
physics<T, SAFE_MATH>::c);
950+
951+
auto w2 = w*w;
952+
auto denome = 1.0 - ec*ec/w2;
953+
auto e11 = 1.0 - (wpe2/w2)/denome;
954+
auto e12 = ((ec/w)*(wpe2/w2))/denome;
955+
auto e33 = wpe2;
956+
957+
for (size_t i = 0, ie = eq->get_num_ion_species(); i < ie; i++) {
958+
const T mi = eq->get_ion_mass(i);
959+
const T charge = static_cast<T> (eq->get_ion_charge(i))
960+
* physics<T, SAFE_MATH>::q;
961+
962+
auto ni = eq->get_ion_density(i, x, y, z);
963+
auto wpi2 = build_plasma_fequency(ni, charge, mi,
964+
physics<T, SAFE_MATH>::c,
965+
physics<T, SAFE_MATH>::epsion0);
966+
auto ic = build_cyclotron_fequency(charge, b_len, mi,
967+
physics<T, SAFE_MATH>::c);
968+
969+
auto denomi = 1.0 - ic*ic/w2;
970+
e11 = e11 - (wpi2/w2)/denomi;
971+
e12 = e12 + ((ic/w)*(wpi2/w2))/denomi;
972+
e33 = e33 + wpi2;
973+
}
974+
e33 = 1.0 - e33/w2;
975+
976+
// Wave numbers.
977+
auto n = (kx*eq->get_esup1(x, y, z) +
978+
ky*eq->get_esup2(x, y, z) +
979+
kz*eq->get_esup3(x, y, z))/w;
980+
auto b_hat = b_vec->unit();
981+
982+
auto npara = b_hat->dot(n);
983+
auto npara2 = npara*npara;
984+
auto nperp = b_hat->cross(n);
985+
auto nperp2 = nperp->dot(nperp);
986+
auto n2 = nperp2 + npara2;
987+
988+
auto es = e11*e11 - e12*e12;
989+
auto A = (e11*nperp2 + e33*npara2)/n2;
990+
auto B = e11*e33*(1.0 + npara2/n2) + es*nperp2/n2;
991+
auto C = e33*es;
992+
993+
const T sign = static_cast<T> (MODE);
994+
return 2.0*A*n2 - B - sign*graph::sqrt(B*B - 4.0*A*C);
995+
}
996+
};
997+
880998
//------------------------------------------------------------------------------
881999
/// @brief Cold Plasma Disperison function.
8821000
///

‎graph_framework/equilibrium.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1360,7 +1360,7 @@ namespace equilibrium {
13601360

13611361
auto q = graph::constant<T, SAFE_MATH> (static_cast<T> (1.60218E-19));
13621362

1363-
ni_cache = te_cache;
1363+
ni_cache = ne_cache;
13641364
ti_cache = (pressure - ne_cache*te_cache*q)/(ni_cache*q);
13651365

13661366
auto phi = graph::atan(x, y);

0 commit comments

Comments
 (0)