Skip to content

Commit 9ad2040

Browse files
cianciosacianciosa
cianciosa
authored and
cianciosa
committed
Add explicit O and X mode dipersion relations.
1 parent c0f6a84 commit 9ad2040

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
@@ -878,6 +878,124 @@ namespace dispersion {
878878
}
879879
};
880880

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

‎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)