Commit 63c3cba6 authored by Cianciosa, Mark's avatar Cianciosa, Mark
Browse files

Add split symplectic integrator to physics tests.

parent b4ca011a

graph_driver/xrays.cpp

+5 −4
Original line number Diff line number Diff line
@@ -39,8 +39,8 @@ int main(int argc, const char * argv[]) { 
    const std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now();



    const size_t num_times = 10000;

    //const size_t num_rays = 1;

    const size_t num_rays = 10000;

    const size_t num_rays = 1;

    //const size_t num_rays = 10000;

    

    std::vector<std::thread> threads(std::max(std::min(std::thread::hardware_concurrency(),

                                                       static_cast<unsigned int> (num_rays)),
@@ -85,8 +85,9 @@ int main(int argc, const char * argv[]) { 
            //solver::split_simplextic<dispersion::bohm_gross<cpu>>

            //solver::rk4<dispersion::bohm_gross<cpu>>

            //solver::rk4<dispersion::simple<cpu>>

            solver::rk4<dispersion::ordinary_wave<cpu>>

            //solver::rk4<dispersion::cold_plasma<cpu>>

            //solver::rk4<dispersion::ordinary_wave<cpu>>

            //solver::rk4<dispersion::extra_ordinary_wave<cpu>>

            solver::rk4<dispersion::cold_plasma<cpu>>

                solve(omega, kx, ky, kz, x, y, z, t, 30.0/num_times, eq);

            solve.init(kx);

            if (thread_number == 0) {

graph_framework.xcodeproj/project.xcworkspace/xcuserdata/m4c.xcuserdatad/UserInterfaceState.xcuserstate

+6.66 KiB (175 KiB)

File changed.

No diff preview for this file type.

View original file
View changed file

graph_framework.xcodeproj/xcuserdata/m4c.xcuserdatad/xcschemes/xcschememanagement.plist

+7 −7
Original line number Diff line number Diff line
@@ -7,17 +7,17 @@ 
		<key>arithmetic_test.xcscheme_^#shared#^_</key>

		<dict>

			<key>orderHint</key>

			<integer>6</integer>

			<integer>8</integer>

		</dict>

		<key>backend_test.xcscheme_^#shared#^_</key>

		<dict>

			<key>orderHint</key>

			<integer>9</integer>

			<integer>7</integer>

		</dict>

		<key>dispersion_test.xcscheme_^#shared#^_</key>

		<dict>

			<key>orderHint</key>

			<integer>4</integer>

			<integer>6</integer>

		</dict>

		<key>graph_driver.xcscheme_^#shared#^_</key>

		<dict>
@@ -27,7 +27,7 @@ 
		<key>graph_framework.xcscheme_^#shared#^_</key>

		<dict>

			<key>orderHint</key>

			<integer>3</integer>

			<integer>5</integer>

		</dict>

		<key>graph_tests.xcscheme_^#shared#^_</key>

		<dict>
@@ -42,7 +42,7 @@ 
		<key>node_test.xcscheme_^#shared#^_</key>

		<dict>

			<key>orderHint</key>

			<integer>7</integer>

			<integer>4</integer>

		</dict>

		<key>physics_test.xcscheme_^#shared#^_</key>

		<dict>
@@ -52,12 +52,12 @@ 
		<key>solver_test.xcscheme_^#shared#^_</key>

		<dict>

			<key>orderHint</key>

			<integer>5</integer>

			<integer>3</integer>

		</dict>

		<key>vector_test.xcscheme_^#shared#^_</key>

		<dict>

			<key>orderHint</key>

			<integer>8</integer>

			<integer>9</integer>

		</dict>

	</dict>

	<key>SuppressBuildableAutocreation</key>

graph_tests/physics_test.cpp

+78 −76
Original line number Diff line number Diff line
@@ -101,50 +101,50 @@ void test_constant() { 
///

///  @param[in] tolarance Tolarance to solver the dispersion function to.

//------------------------------------------------------------------------------

template<typename BACKEND>

void test_bohm_gross(const typename BACKEND::base tolarance) {

template<typename SOLVER>

void test_bohm_gross(const typename SOLVER::base tolarance) {

    std::mt19937_64 engine(static_cast<uint64_t> (std::chrono::system_clock::to_time_t(std::chrono::system_clock::now())));

    std::uniform_real_distribution<double> real_dist(0.1, 1.0);



    auto omega = graph::variable<BACKEND> (1, "\\omega");

    auto kx = graph::variable<BACKEND> (1, "k_{x}");

    auto ky = graph::variable<BACKEND> (1, "k_{y}");

    auto kz = graph::variable<BACKEND> (1, "k_{z}");

    auto x = graph::variable<BACKEND> (1, "x");

    auto y = graph::variable<BACKEND> (1, "y");

    auto z = graph::variable<BACKEND> (1, "z");

    auto t = graph::variable<BACKEND> (1, "t");

    auto omega = graph::variable<typename SOLVER::backend> (1, "\\omega");

    auto kx = graph::variable<typename SOLVER::backend> (1, "k_{x}");

    auto ky = graph::variable<typename SOLVER::backend> (1, "k_{y}");

    auto kz = graph::variable<typename SOLVER::backend> (1, "k_{z}");

    auto x = graph::variable<typename SOLVER::backend> (1, "x");

    auto y = graph::variable<typename SOLVER::backend> (1, "y");

    auto z = graph::variable<typename SOLVER::backend> (1, "z");

    auto t = graph::variable<typename SOLVER::backend> (1, "t");



//  Constants

    const typename BACKEND::base q = 1.602176634E-19;

    const typename BACKEND::base me = 9.1093837015E-31;

    const typename BACKEND::base mu0 = M_PI*4.0E-7;

    const typename BACKEND::base epsilon0 = 8.8541878138E-12;

    const typename BACKEND::base c = 1.0/sqrt(mu0*epsilon0);

    const typename BACKEND::base omega0 = 600.0;

    const typename BACKEND::base ne0 = 1.0E19;

    const typename BACKEND::base te = 1000.0;

    const typename SOLVER::base q = 1.602176634E-19;

    const typename SOLVER::base me = 9.1093837015E-31;

    const typename SOLVER::base mu0 = M_PI*4.0E-7;

    const typename SOLVER::base epsilon0 = 8.8541878138E-12;

    const typename SOLVER::base c = 1.0/sqrt(mu0*epsilon0);

    const typename SOLVER::base omega0 = 600.0;

    const typename SOLVER::base ne0 = 1.0E19;

    const typename SOLVER::base te = 1000.0;

    

    const typename BACKEND::base omega2 = (ne0*0.9*q*q)/(epsilon0*me*c*c);

    const typename BACKEND::base omega2p = (ne0*0.1*q*q)/(epsilon0*me*c*c);

    const typename BACKEND::base vth2 = 2*1.602176634E-19*te/(me*c*c);

    const typename SOLVER::base omega2 = (ne0*0.9*q*q)/(epsilon0*me*c*c);

    const typename SOLVER::base omega2p = (ne0*0.1*q*q)/(epsilon0*me*c*c);

    const typename SOLVER::base vth2 = 2*1.602176634E-19*te/(me*c*c);

    

    const typename BACKEND::base k0 = std::sqrt(2.0/3.0*(omega0*omega0 - omega2)/vth2);

    const typename SOLVER::base k0 = std::sqrt(2.0/3.0*(omega0*omega0 - omega2)/vth2);

    

//  Omega must be greater than plasma frequency for the wave to propagate.

    omega->set(backend::base_cast<BACKEND> (600.0));

    kx->set(backend::base_cast<BACKEND> (1000.0));

    ky->set(backend::base_cast<BACKEND> (0.0));

    kz->set(backend::base_cast<BACKEND> (0.0));

    x->set(backend::base_cast<BACKEND> (-1.0));

    y->set(backend::base_cast<BACKEND> (0.0));

    z->set(backend::base_cast<BACKEND> (0.0));

    t->set(backend::base_cast<BACKEND> (0.0));



    const typename BACKEND::base dt = 0.1;

    

    auto eq = equilibrium::make_no_magnetic_field<BACKEND> ();

    solver::rk4<dispersion::bohm_gross<BACKEND>> solve(omega, kx, ky, kz, x, y, z, t, dt, eq);

    omega->set(backend::base_cast<typename SOLVER::backend> (600.0));

    kx->set(backend::base_cast<typename SOLVER::backend> (1000.0));

    ky->set(backend::base_cast<typename SOLVER::backend> (0.0));

    kz->set(backend::base_cast<typename SOLVER::backend> (0.0));

    x->set(backend::base_cast<typename SOLVER::backend> (-1.0));

    y->set(backend::base_cast<typename SOLVER::backend> (0.0));

    z->set(backend::base_cast<typename SOLVER::backend> (0.0));

    t->set(backend::base_cast<typename SOLVER::backend> (0.0));



    const typename SOLVER::base dt = 0.1;

    

    auto eq = equilibrium::make_no_magnetic_field<typename SOLVER::backend> ();

    SOLVER solve(omega, kx, ky, kz, x, y, z, t, dt, eq);

    solve.init(kx);

    

    const auto diff = kx->evaluate().at(0) - k0;
@@ -154,8 +154,8 @@ void test_bohm_gross(const typename BACKEND::base tolarance) { 
    for (size_t i = 0; i < 20; i++) {

        solve.step();

    }

    const typename BACKEND::base time = t->evaluate().at(0);

    const typename BACKEND::base expected_x = -3.0/8.0*vth2*omega2p/(omega0*omega0)*time*time

    const typename SOLVER::base time = t->evaluate().at(0);

    const typename SOLVER::base expected_x = -3.0/8.0*vth2*omega2p/(omega0*omega0)*time*time

                                           + 3.0/2.0*vth2/omega0*k0*time - 1.0;

    

    const auto diff_x = x->evaluate().at(0) - expected_x;
@@ -198,48 +198,48 @@ void test_bohm_gross(const typename BACKEND::base tolarance) { 
///

///  @param[in] tolarance Tolarance to solver the dispersion function to.

//------------------------------------------------------------------------------

template<typename BACKEND>

void test_light_wave(const typename BACKEND::base tolarance) {

template<typename SOLVER>

void test_light_wave(const typename SOLVER::base tolarance) {

    std::mt19937_64 engine(static_cast<uint64_t> (std::chrono::system_clock::to_time_t(std::chrono::system_clock::now())));

    std::uniform_real_distribution<double> real_dist(0.1, 1.0);



    auto omega = graph::variable<BACKEND> (1, "\\omega");

    auto kx = graph::variable<BACKEND> (1, "k_{x}");

    auto ky = graph::variable<BACKEND> (1, "k_{y}");

    auto kz = graph::variable<BACKEND> (1, "k_{z}");

    auto x = graph::variable<BACKEND> (1, "x");

    auto y = graph::variable<BACKEND> (1, "y");

    auto z = graph::variable<BACKEND> (1, "z");

    auto t = graph::variable<BACKEND> (1, "t");

    auto omega = graph::variable<typename SOLVER::backend> (1, "\\omega");

    auto kx = graph::variable<typename SOLVER::backend> (1, "k_{x}");

    auto ky = graph::variable<typename SOLVER::backend> (1, "k_{y}");

    auto kz = graph::variable<typename SOLVER::backend> (1, "k_{z}");

    auto x = graph::variable<typename SOLVER::backend> (1, "x");

    auto y = graph::variable<typename SOLVER::backend> (1, "y");

    auto z = graph::variable<typename SOLVER::backend> (1, "z");

    auto t = graph::variable<typename SOLVER::backend> (1, "t");



//  Constants

    const typename BACKEND::base q = 1.602176634E-19;

    const typename BACKEND::base me = 9.1093837015E-31;

    const typename BACKEND::base mu0 = M_PI*4.0E-7;

    const typename BACKEND::base epsilon0 = 8.8541878138E-12;

    const typename BACKEND::base c = 1.0/sqrt(mu0*epsilon0);

    const typename BACKEND::base omega0 = 600.0;

    const typename BACKEND::base ne0 = 1.0E19;

    const typename SOLVER::base q = 1.602176634E-19;

    const typename SOLVER::base me = 9.1093837015E-31;

    const typename SOLVER::base mu0 = M_PI*4.0E-7;

    const typename SOLVER::base epsilon0 = 8.8541878138E-12;

    const typename SOLVER::base c = 1.0/sqrt(mu0*epsilon0);

    const typename SOLVER::base omega0 = 600.0;

    const typename SOLVER::base ne0 = 1.0E19;

    

    const typename BACKEND::base omega2 = (ne0*0.9*q*q)/(epsilon0*me*c*c);

    const typename BACKEND::base omega2p = (ne0*0.1*q*q)/(epsilon0*me*c*c);

    const typename SOLVER::base omega2 = (ne0*0.9*q*q)/(epsilon0*me*c*c);

    const typename SOLVER::base omega2p = (ne0*0.1*q*q)/(epsilon0*me*c*c);

    

    const typename BACKEND::base k0 = std::sqrt(omega0*omega0 - omega2);

    const typename SOLVER::base k0 = std::sqrt(omega0*omega0 - omega2);

    

//  Omega must be greater than plasma frequency for the wave to propagate.

    omega->set(backend::base_cast<BACKEND> (600.0));

    kx->set(backend::base_cast<BACKEND> (100.0));

    ky->set(backend::base_cast<BACKEND> (0.0));

    kz->set(backend::base_cast<BACKEND> (0.0));

    x->set(backend::base_cast<BACKEND> (-1.0));

    y->set(backend::base_cast<BACKEND> (0.0));

    z->set(backend::base_cast<BACKEND> (0.0));

    t->set(backend::base_cast<BACKEND> (0.0));



    const typename BACKEND::base dt = 0.1;

    

    auto eq = equilibrium::make_no_magnetic_field<BACKEND> ();

    solver::rk4<dispersion::light_wave<BACKEND>> solve(omega, kx, ky, kz, x, y, z, t, dt, eq);

    omega->set(backend::base_cast<typename SOLVER::backend> (600.0));

    kx->set(backend::base_cast<typename SOLVER::backend> (100.0));

    ky->set(backend::base_cast<typename SOLVER::backend> (0.0));

    kz->set(backend::base_cast<typename SOLVER::backend> (0.0));

    x->set(backend::base_cast<typename SOLVER::backend> (-1.0));

    y->set(backend::base_cast<typename SOLVER::backend> (0.0));

    z->set(backend::base_cast<typename SOLVER::backend> (0.0));

    t->set(backend::base_cast<typename SOLVER::backend> (0.0));



    const typename SOLVER::base dt = 0.1;

    

    auto eq = equilibrium::make_no_magnetic_field<typename SOLVER::backend> ();

    SOLVER solve(omega, kx, ky, kz, x, y, z, t, dt, eq);

    solve.init(kx, tolarance);

    

    const auto diff = kx->evaluate().at(0) - k0;
@@ -249,8 +249,8 @@ void test_light_wave(const typename BACKEND::base tolarance) { 
    for (size_t i = 0; i < 20; i++) {

        solve.step();

    }

    const typename BACKEND::base time = t->evaluate().at(0);

    const typename BACKEND::base expected_x = -omega2p/(4.0*omega0*omega0)*time*time

    const typename SOLVER::base time = t->evaluate().at(0);

    const typename SOLVER::base expected_x = -omega2p/(4.0*omega0*omega0)*time*time

                                           + k0/omega0*time - 1.0;

    

    const auto diff_x = x->evaluate().at(0) - expected_x;
@@ -570,8 +570,10 @@ void test_reflection(const typename BACKEND::base tolarance, 
//------------------------------------------------------------------------------

template<typename BACKEND> void run_tests(const typename BACKEND::base tolarance) {

    test_constant<BACKEND> ();

    test_bohm_gross<BACKEND> (tolarance);

    test_light_wave<BACKEND> (tolarance);

    test_bohm_gross<solver::rk4<dispersion::bohm_gross<BACKEND>>> (tolarance);

    test_bohm_gross<solver::split_simplextic<dispersion::bohm_gross<BACKEND>>> (tolarance);

    test_light_wave<solver::rk4<dispersion::light_wave<BACKEND>>> (tolarance);

    test_light_wave<solver::split_simplextic<dispersion::light_wave<BACKEND>>> (tolarance);

    test_acoustic_wave<BACKEND> (tolarance);

    test_o_mode_wave<BACKEND> ();

    test_reflection<BACKEND> (tolarance, 0.7, 0.1, 22.0);

graph_framework/solver.hpp

+9 −9

File changed.

Contains only whitespace changes.