Commit d33632ca authored by Cianciosa, Mark's avatar Cianciosa, Mark
Browse files

Fix sign of electron cyclotron frequency.

parent cd5cfe90
Loading
Loading
Loading
Loading
+18 −17
Original line number Diff line number Diff line
@@ -11,8 +11,8 @@
#include "../graph_framework/solver.hpp"
#include "../graph_framework/timing.hpp"

const bool print = false;
const bool write_step = true;
const bool print = true;
const bool write_step = false;
const bool print_expressions = false;

//------------------------------------------------------------------------------
@@ -26,19 +26,19 @@ int main(int argc, const char * argv[]) {

    std::mutex sync;

    typedef float base;
    //typedef float base;
    //typedef double base;
    //typedef std::complex<float> base;
    //typedef std::complex<double> base;
    //constexpr bool use_safe_math = true;
    constexpr bool use_safe_math = false;
    typedef std::complex<double> base;
    constexpr bool use_safe_math = true;
    //constexpr bool use_safe_math = false;

    const timeing::measure_diagnostic total("Total Time");

    const size_t num_times = 10000;
    const size_t sub_steps = 1;
    const size_t num_times = 100000;
    const size_t sub_steps = 10;
    const size_t num_steps = num_times/sub_steps;
    const size_t num_rays = 1000000;
    const size_t num_rays = 1;

    std::vector<std::thread> threads(0);
    if constexpr (jit::use_gpu<base> ()) {
@@ -82,21 +82,22 @@ int main(int argc, const char * argv[]) {
                }
            }

            x->set(static_cast<base> (2.5));
            omega->set(static_cast<base> (500.0));
            x->set(static_cast<base> (-12.0));
            //x->set(static_cast<base> (0.0));
            y->set(static_cast<base> (0.0));
            z->set(static_cast<base> (0.0));
            kx->set(static_cast<base> (-600.0));
            kx->set(static_cast<base> (-600));
            //kx->set(static_cast<base> (600.0));
            ky->set(static_cast<base> (0.0));
            kz->set(static_cast<base> (0.0));
            kz->set(static_cast<base> (10.0));

            auto eq = equilibrium::make_efit<base, use_safe_math> (NC_FILE, sync);
            //auto eq = equilibrium::make_efit<base, use_safe_math> (NC_FILE, sync);
            //auto eq = equilibrium::make_slab_density<base, use_safe_math> ();
            //auto eq = equilibrium::make_slab_field<base, use_safe_math> ();
            auto eq = equilibrium::make_slab_field<base, use_safe_math> ();
            //auto eq = equilibrium::make_no_magnetic_field<base, use_safe_math> ();

            const base endtime = static_cast<base> (1.0);
            const base endtime = static_cast<base> (10.0);
            //const base endtime = static_cast<base> (0.25);
            const base dt = endtime/static_cast<base> (num_times);

@@ -111,9 +112,9 @@ int main(int argc, const char * argv[]) {
            //solver::rk4<dispersion::simple<base, use_safe_math>>
            //solver::rk4<dispersion::ordinary_wave<base, use_safe_math>>
            //solver::rk4<dispersion::extra_ordinary_wave<base, use_safe_math>>
            solver::rk4<dispersion::cold_plasma<base, use_safe_math>>
            //solver::rk4<dispersion::cold_plasma<base, use_safe_math>>
            //solver::adaptive_rk4<dispersion::ordinary_wave<base, use_safe_math>>
            //solver::rk4<dispersion::hot_plasma<base, dispersion::z_erfi<base, use_safe_math>, use_safe_math>>
            solver::rk4<dispersion::hot_plasma<base, dispersion::z_erfi<base, use_safe_math>, use_safe_math>>
            //solver::rk4<dispersion::hot_plasma_expandion<base, dispersion::z_erfi<base, use_safe_math>, use_safe_math>>
                solve(omega, kx, ky, kz, x, y, z, t, dt, eq,
                      stream.str(), local_num_rays);
+2 −2
Original line number Diff line number Diff line
@@ -1077,7 +1077,7 @@ namespace dispersion {
                                  physics<T, SAFE_MATH>::me);

//  Setup characteristic frequencies.
            auto ec = build_cyclotron_fequency(none*physics<T, SAFE_MATH>::q,
            auto ec = build_cyclotron_fequency(physics<T, SAFE_MATH>::q,
                                               b_len,
                                               physics<T, SAFE_MATH>::me,
                                               physics<T, SAFE_MATH>::c);
@@ -1198,7 +1198,7 @@ namespace dispersion {
            auto ve = graph::sqrt(two*physics<T, SAFE_MATH>::q*te/physics<T, SAFE_MATH>::me);

//  Setup characteristic frequencies.
            auto ec = build_cyclotron_fequency(none*physics<T, SAFE_MATH>::q, b_len,
            auto ec = build_cyclotron_fequency(physics<T, SAFE_MATH>::q, b_len,
                                               physics<T, SAFE_MATH>::me,
                                               physics<T, SAFE_MATH>::c);
            auto wpe2 = build_plasma_fequency(ne, physics<T, SAFE_MATH>::q,
+4 −4
Original line number Diff line number Diff line
@@ -487,7 +487,7 @@ namespace equilibrium {
                             graph::shared_leaf<T, SAFE_MATH> y,
                             graph::shared_leaf<T, SAFE_MATH> z) final {
            return graph::constant<T, SAFE_MATH> (static_cast<T> (1.0E19)) *
                   (graph::constant<T, SAFE_MATH> (static_cast<T> (0.1))*x +
                   (graph::constant<T, SAFE_MATH> (static_cast<T> (0.01))*x +
                    graph::one<T, SAFE_MATH> ());
        }

@@ -517,8 +517,8 @@ namespace equilibrium {
        get_electron_temperature(graph::shared_leaf<T, SAFE_MATH> x,
                                 graph::shared_leaf<T, SAFE_MATH> y,
                                 graph::shared_leaf<T, SAFE_MATH> z) final {
            return graph::constant<T, SAFE_MATH> (static_cast<T> (1000.0)) *
                   (graph::constant<T, SAFE_MATH> (static_cast<T> (0.1))*x +
            return graph::constant<T, SAFE_MATH> (static_cast<T> (2000.0)) *
                   (graph::constant<T, SAFE_MATH> (static_cast<T> (0.01))*x +
                    graph::one<T, SAFE_MATH> ());
        }

@@ -550,7 +550,7 @@ namespace equilibrium {
                           graph::shared_leaf<T, SAFE_MATH> z) final {
            auto zero = graph::zero<T, SAFE_MATH> ();
            return graph::vector(zero, zero,
                                 graph::constant<T, SAFE_MATH> (static_cast<T> (0.1))*x +
                                 graph::constant<T, SAFE_MATH> (static_cast<T> (0.01))*x +
                                 graph::one<T, SAFE_MATH> ());
        }
    };