Merge pull request #161 from taiga-project/feature-157

Feature 157
taiga-project · May 9, 2022 · 5848283 · 5848283
2 parents fc5c75b + ee4b175
commit 5848283
Show file tree

Hide file tree

Showing 7 changed files with 356 additions and 5 deletions.
diff --git a/Makefile b/Makefile
@@ -47,6 +47,12 @@ test_init: tests | $(BIN)
 test_framework: tests  | $(BIN)
 	$(GCC) $(DEFAULT_FLAGS) -o $(BIN)/test_framework.exe tests/taiga_test_example.c
 
+example_solvers: $(OBJ)/example_solvers.o $(OBJ)/test_solver.o | $(BIN)
+	$(GCC) $(DEFAULT_FLAGS) $^ -lm -o $(BIN)/example_solvers.exe
+
+$(OBJ)/example_solvers.o: example/solvers/export_solver.c $(OBJ)
+	$(GCC) $(DEFAULT_FLAGS) -w -Isrc -Iinclude -Itests -c $< -lm -o $@
+
 field: tests | $(BIN)
 	$(NVCC) $(CFLAGS) $(DEFAULT_FLAGS) -o $(BIN)/test_field.exe tests/test_field.cu
 

diff --git a/example/solvers/__init__.py b/example/solvers/__init__.py
diff --git a/example/solvers/compare_solvers.py b/example/solvers/compare_solvers.py
@@ -0,0 +1,140 @@
+import numpy
+import matplotlib.pyplot
+
+
+class SolverData:
+    directory = "bin/example"
+    reference_velocity = 400000
+    number_of_steps = 20000000
+    frequency_of_export = 100000
+    timestep = 1e-9
+    field_scenario = 'B = 1'
+
+    def __init__(self, solver_id, scenario='default'):
+        self.solver_file = solver_id
+        self.solver_name = self.get_solver_name(solver_id)
+        self.solver_file += '_' + scenario
+        self.set_scenario(scenario)
+        self.data = self.get_data()
+        self.x = self.get_field(0)
+        self.y = self.get_field(1)
+        self.v = numpy.sqrt(self.get_field(3)**2 + self.get_field(4)**2 + self.get_field(5)**2)
+        self.v_relative = self.v / self.reference_velocity
+        self.x_axis = numpy.arange(0, self.number_of_steps, self.frequency_of_export)
+        self.line_style = self.get_line_style(solver_id)
+
+    def set_scenario(self, scenario):
+        if scenario == 'start10':
+            self.timestep = 1e-10
+
+        if scenario in {'start8', 'gradb_start8', 'b__r_start8'}:
+            self.timestep = 1e-8
+
+        if scenario in {'gradb', 'gradb_start'}:
+            self.field_scenario = 'B_z = 1 + 0.01 y'
+
+        if scenario in {'eparb', 'eparb_start'}:
+            self.field_scenario = 'B_z = 1, E_z = 0.01'
+
+        if scenario in {'b__r', 'b__r_start', 'b__r_start8'}:
+            self.field_scenario = '$B_z = \sqrt{x^2+y^2}$'
+
+        if scenario in {'start', 'gradb_start', 'b__r_start', 'eparb_start'}:
+            self.number_of_steps = 100000
+            self.frequency_of_export = 10
+
+        if scenario in {'start8', 'gradb_start8', 'b__r_start8'}:
+            self.number_of_steps = 10000
+            self.frequency_of_export = 1
+
+    @staticmethod
+    def get_solver_name(solver_id):
+        if solver_id == 'rk4':
+            return 'linearised Runge--Kutta'
+        elif solver_id == 'rkn':
+            return 'Runge--Kutta--Nyström'
+        elif solver_id == 'verlet':
+            return 'Verlet--Boris'
+        elif solver_id == 'yoshida':
+            return 'Yoshida--Boris'
+        else:
+            return 'unknown'
+
+    @staticmethod
+    def get_line_style(solver_id):
+        if solver_id == 'rk4':
+            return '-'
+        elif solver_id == 'rkn':
+            return '--'
+        elif solver_id == 'verlet':
+            return '-'
+        elif solver_id == 'yoshida':
+            return '--'
+        else:
+            return ':'
+
+    def get_data(self):
+        return numpy.loadtxt(self.directory + '/' + self.solver_file + '.dat')
+
+    def get_field(self, identifier):
+        return self.data[:, identifier]
+
+    def get_timestep(self):
+        exponent = numpy.floor(numpy.log10(self.timestep))
+        mantissa = self.timestep / 10 ** exponent
+        return '$' + (str(mantissa) + ' \cdot ' if (mantissa != 1.0) else '') + '10^{' + str(int(exponent)) + '}$ s'
+
+    def plot_xy(self):
+        matplotlib.pyplot.plot(self.x, self.y, '.', markersize=1.5, label=self.solver_name)
+        matplotlib.pyplot.title('timestep: ' + self.get_timestep() + ', ' + self.field_scenario)
+
+    def plot_v(self):
+        matplotlib.pyplot.plot(self.x_axis, self.v_relative, self.line_style, linewidth=2.5, label=self.solver_name)
+        matplotlib.pyplot.title('timestep: ' + self.get_timestep() + ', ' + self.field_scenario)
+
+
+def visualise():
+    matplotlib.pyplot.xlabel('steps')
+    matplotlib.pyplot.ylabel('relative velocity')
+    matplotlib.pyplot.legend()
+    matplotlib.pyplot.tight_layout()
+
+
+def visualise_trajectory():
+    matplotlib.pyplot.xlabel('x [m]')
+    matplotlib.pyplot.ylabel('y [m]')
+    matplotlib.pyplot.axis('square')
+    matplotlib.pyplot.legend()
+    matplotlib.pyplot.tight_layout()
+
+
+def compare_solvers(scenario='default'):
+    SolverData('rk4', scenario).plot_v()
+    SolverData('rkn', scenario).plot_v()
+    SolverData('verlet', scenario).plot_v()
+    SolverData('yoshida', scenario).plot_v()
+    visualise()
+    matplotlib.pyplot.savefig('bin/example/' + scenario + '.pdf')
+    matplotlib.pyplot.savefig('bin/example/' + scenario + '.png')
+    matplotlib.pyplot.show()
+
+
+def compare_solvers_trajectory(scenario='default'):
+    SolverData('rk4', scenario).plot_xy()
+    SolverData('rkn', scenario).plot_xy()
+    SolverData('verlet', scenario).plot_xy()
+    SolverData('yoshida', scenario).plot_xy()
+    visualise_trajectory()
+    matplotlib.pyplot.savefig('bin/example/' + scenario + '_traj.pdf')
+    matplotlib.pyplot.savefig('bin/example/' + scenario + '_traj.png')
+    matplotlib.pyplot.show()
+
+
+if __name__ == "__main__":
+    compare_solvers('b__r')
+    compare_solvers_trajectory('b__r')
+    compare_solvers('start8')
+    compare_solvers_trajectory('start8')
+    compare_solvers_trajectory('start')
+    compare_solvers('gradb_start')
+    compare_solvers_trajectory('gradb_start')
diff --git a/example/solvers/export_solver.c b/example/solvers/export_solver.c
@@ -0,0 +1,153 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdbool.h>
+#include <math.h>
+#include <string.h>
+#include <sys/stat.h>
+
+#include "export_solver.h"
+#include "test_solver.h"
+
+#include "utils/taiga_constants.h"
+#include "utils/prop.h"
+#include "core/maths/maths.cuh"
+#include "core/solvers/solvers.cuh"
+#include "core/solvers/rk4.cuh"
+#include "core/solvers/runge_kutta_nystrom.cuh"
+#include "core/solvers/yoshida.cuh"
+#include "core/solvers/verlet.cuh"
+#include "core/solvers/boris.cuh"
+#include "core/physics/lorentz.cuh"
+
+#include "utils/basic_functions.c"
+
+#define FOLDER "example"
+
+void export_coordinate (FILE *f, double *X, long N_half) {
+    fprintf(f, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%ld\n", X[0], X[1], X[2], X[3], X[4], X[5], N_half);
+}
+
+void generate_R_B_const_field(double *X, double *local_bfield, double *local_efield,
+                              TaigaCommons *c, bool is_electric_field_on,
+                              int *local_spline_indices,
+                              double *local_spline_brad, double *local_spline_bz, double *local_spline_btor,
+                              double *local_spline_erad, double *local_spline_ez, double *local_spline_etor,
+                              double *local_psi_n) {
+    local_bfield[0] = 0.0;
+    local_bfield[1] = 0.0;
+    local_bfield[2] = 1.0 / (X[0] + 1.0);
+    local_efield[0] = 0.0;
+    local_efield[1] = 0.0;
+    local_efield[2] = 0.0;
+}
+
+void generate_B_over_R_field(double *X, double *local_bfield, double *local_efield,
+                             TaigaCommons *c, bool is_electric_field_on,
+                             int *local_spline_indices,
+                             double *local_spline_brad, double *local_spline_bz, double *local_spline_btor,
+                             double *local_spline_erad, double *local_spline_ez, double *local_spline_etor,
+                             double *local_psi_n) {
+    local_bfield[0] = 0.0;
+    local_bfield[1] = 0.0;
+    local_bfield[2] = sqrt (X[0] * X[0] + X[1] * X[1]);
+    local_efield[0] = 0.0;
+    local_efield[1] = 0.0;
+    local_efield[2] = 0.0;
+}
+
+void run_field_with_solver_and_export(char* scenario_name, double timestep, int field_type, char* solver_name,
+                                      long number_of_steps, long frequency_of_export,
+                                      double (*solve_diffeq)(double *X, double eperm, double timestep,
+                                                             TaigaCommons *c, bool is_electric_field_on,
+                                                             int *local_spline_indices,
+                                                             double *local_spline_brad, double *local_spline_bz, double *local_spline_btor,
+                                                             double *local_spline_erad, double *local_spline_ez, double *local_spline_etor,
+                                                             double *local_psi_n) ) {
+    switch(field_type){
+        case HOMOGENEOUS:
+            generate_local_field = &generate_homogeneous_magnetic_field;
+            break;
+        case GRAD_B:
+            generate_local_field = &generate_grad_B_field;
+            break;
+        case E_PAR_B:
+            generate_local_field = &generate_E_par_B_field;
+            break;
+        case BR_FIELD:
+            generate_local_field = &generate_R_B_const_field;
+            break;
+        case B_OVER_R_FIELD:
+            generate_local_field = &generate_B_over_R_field;
+            break;
+        default:
+            printf("Error: Illegal field_type\n");
+    }
+    FILE *file;
+    double X[6] = {0};
+    int local_spline_indices[2];
+    local_spline_indices[0] = SPLINE_INDEX_ERROR;
+    local_spline_indices[1] = SPLINE_INDEX_ERROR;
+    double local_spline_brad[16];
+    double local_spline_bz[16];
+    double local_spline_btor[16];
+    double local_spline_erad[16];
+    double local_spline_ez[16];
+    double local_spline_etor[16];
+    double local_psi_n[16];
+    TaigaCommons *c;
+    double eperm = 4e7;
+    int is_electric_field_on = true;
+    long number_of_half_cyclotron_periods = 0;
+    double v_previous;
+
+    get_acceleration_from_lorentz_force = &get_acceleration_from_lorentz_force_with_electric_field;
+    X[4] = eperm * LARMOR_RADIUS;
+
+    if (field_type == B_OVER_R_FIELD)   X[0] = 1.5;
+
+    mkdir(FOLDER, S_IRWXU | S_IRWXG | S_IRWXO);
+    const char *path = concat(FOLDER, "/", solver_name, "_", scenario_name, ".dat", NULL);
+    printf("export to: %s\n", path);
+    file = fopen(path ,"w");
+    for (int i = 0; i < number_of_steps; ++i) {
+        v_previous = X[4];
+        solve_diffeq(X, eperm, timestep,
+                     c, is_electric_field_on,
+                     local_spline_indices,
+                     local_spline_brad, local_spline_bz, local_spline_btor,
+                     local_spline_erad, local_spline_ez, local_spline_etor,
+                     local_psi_n);
+        if ((X[4] * v_previous) <= 0)   ++number_of_half_cyclotron_periods;
+        if (i % frequency_of_export == 0)   export_coordinate(file, X, number_of_half_cyclotron_periods);
+    }
+    fclose(file);
+}
+
+
+void run_scenario(char* scenario_name, double timestep, int field_type,
+                  long number_of_steps, long frequency_of_export){
+    run_field_with_solver_and_export(scenario_name, timestep, field_type, "rk4", number_of_steps, frequency_of_export, solve_diffeq_by_rk4);
+    run_field_with_solver_and_export(scenario_name, timestep, field_type, "rkn", number_of_steps, frequency_of_export, solve_diffeq_by_rkn);
+    run_field_with_solver_and_export(scenario_name, timestep, field_type, "verlet", number_of_steps, frequency_of_export, solve_diffeq_by_verlet);
+    run_field_with_solver_and_export(scenario_name, timestep, field_type, "yoshida", number_of_steps, frequency_of_export, solve_diffeq_by_yoshida);
+}
+
+int main() {
+    double timestep = 1e-9;
+    long number_of_steps = 20000000;
+    long frequency_of_export = 100000;
+
+    run_scenario("default", timestep, HOMOGENEOUS, number_of_steps, frequency_of_export);
+    run_scenario("gradb", timestep, GRAD_B, number_of_steps, frequency_of_export);
+    run_scenario("eparb", timestep, E_PAR_B, number_of_steps, frequency_of_export);
+    run_scenario("br", timestep, BR_FIELD, number_of_steps, frequency_of_export);
+    run_scenario("b__r", timestep, B_OVER_R_FIELD, number_of_steps, frequency_of_export);
+    run_scenario("start8", 1e-8, HOMOGENEOUS, 10000, 1);
+    run_scenario("b__r_start8", 1e-8, B_OVER_R_FIELD, 10000, 1);
+    run_scenario("start", timestep, HOMOGENEOUS, 100000, 10);
+    run_scenario("b__r_start", timestep, B_OVER_R_FIELD, 100000, 10);
+    run_scenario("gradb_start", timestep, GRAD_B, 100000, 10);
+    run_scenario("eparb_start", timestep, E_PAR_B, 100000, 10);
+
+    return 0;
+}
diff --git a/example/solvers/export_solver.h b/example/solvers/export_solver.h
@@ -0,0 +1,37 @@
+#ifndef EXPORT_SOLVER_CUH
+#define EXPORT_SOLVER_CUH
+
+#include <stdbool.h>
+#include "utils/taiga_constants.h"
+#include "utils/prop.h"
+
+#define __device__ ;
+
+#define BR_FIELD 90
+#define B_OVER_R_FIELD 91
+
+void export_coordinate (FILE *f, double *X, long N_half);
+void generate_R_B_const_field(double *X, double *local_bfield, double *local_efield,
+                              TaigaCommons *c, bool is_electric_field_on,
+                              int *local_spline_indices,
+                              double *local_spline_brad, double *local_spline_bz, double *local_spline_btor,
+                              double *local_spline_erad, double *local_spline_ez, double *local_spline_etor,
+                              double *local_psi_n);
+void generate_B_over_R_field(double *X, double *local_bfield, double *local_efield,
+                             TaigaCommons *c, bool is_electric_field_on,
+                             int *local_spline_indices,
+                             double *local_spline_brad, double *local_spline_bz, double *local_spline_btor,
+                             double *local_spline_erad, double *local_spline_ez, double *local_spline_etor,
+                             double *local_psi_n);
+void run_field_with_solver_and_export(char* scenario_name, double timestep, int field_type, char* solver_name,
+                                      long number_of_steps, long frequency_of_export,
+                                      double (*solve_diffeq)(double *X, double eperm, double timestep,
+                                                             TaigaCommons *c, bool is_electric_field_on,
+                                                             int *local_spline_indices,
+                                                             double *local_spline_brad, double *local_spline_bz, double *local_spline_btor,
+                                                             double *local_spline_erad, double *local_spline_ez, double *local_spline_etor,
+                                                             double *local_psi_n) );
+void run_scenario(char* scenario_name, double timestep, int field_type,
+                  long number_of_steps, long frequency_of_export);
+
+#endif //EXPORT_SOLVER_CUH
diff --git a/example/solvers/test_compare_solvers.py b/example/solvers/test_compare_solvers.py
@@ -0,0 +1,15 @@
+import unittest
+from taiga.example.solvers.compare_solvers import SolverData
+
+
+class TestCompareSolvers(unittest):
+    def test_timestep(self):
+        s = SolverData('rk4')
+        reference = '$10^{-9}$ s'
+        self.assertEqual(reference, s.get_timestep())
+
+    def test_timestep_1p5em10(self):
+        s = SolverData('rk4')
+        s.timestep = 1.5e-10
+        reference = '$1.5 \cdot 10^{-10}$ s'
+        self.assertEqual(reference, s.get_timestep())
diff --git a/tests/test_solver.c b/tests/test_solver.c
@@ -18,7 +18,7 @@
 #include "core/solvers/boris.cu"
 #include "core/physics/lorentz.cu"
 
-#define GRAD_B_FACTOR 0.0001
+#define GRAD_B_FACTOR 0.01
 #define E_OVER_B 1
 
 #define NUMBER_OF_CYCLOTRON_PERIODS 1000
@@ -189,10 +189,10 @@ int test_solver() {
 
     double reference_grad_B_drift = -GRAD_B_FACTOR * PI * LARMOR_RADIUS * LARMOR_RADIUS * NUMBER_OF_CYCLOTRON_PERIODS_GRAD_B;
     double reference_speed = eperm * LARMOR_RADIUS;
-    TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_grad_B_drift, run_field_with_solver(timestep, GRAD_B, GET_POSITION, solve_diffeq_by_rk4), 1e-5, "4th order linearised Runge--Kutta (grad B)");
-    TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_grad_B_drift, run_field_with_solver(timestep, GRAD_B, GET_POSITION, solve_diffeq_by_rkn), 1e-5, "4th order linearised Runge--Kutta--Nystrom (grad B)");
-    TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_grad_B_drift, run_field_with_solver(timestep, GRAD_B, GET_POSITION, solve_diffeq_by_verlet), 1e-5, "velocity-Verlet--Boris (grad B)");
-    TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_grad_B_drift, run_field_with_solver(timestep, GRAD_B, GET_POSITION, solve_diffeq_by_yoshida), 1e-5, "Yoshida--Boris (grad B)");
+    TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_grad_B_drift, run_field_with_solver(timestep, GRAD_B, GET_POSITION, solve_diffeq_by_rk4), 3e-5, "4th order linearised Runge--Kutta (grad B)");
+    TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_grad_B_drift, run_field_with_solver(timestep, GRAD_B, GET_POSITION, solve_diffeq_by_rkn), 3e-5, "4th order linearised Runge--Kutta--Nystrom (grad B)");
+    TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_grad_B_drift, run_field_with_solver(timestep, GRAD_B, GET_POSITION, solve_diffeq_by_verlet), 3e-5, "velocity-Verlet--Boris (grad B)");
+    TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_grad_B_drift, run_field_with_solver(timestep, GRAD_B, GET_POSITION, solve_diffeq_by_yoshida), 3e-5, "Yoshida--Boris (grad B)");
 
     TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_speed, run_field_with_solver(timestep, GRAD_B, GET_SPEED, solve_diffeq_by_rk4), 1000, "4th order linearised Runge--Kutta (grad B, speed)");
     TAIGA_ASSERT_ALMOST_EQ_MAX_DIFF(reference_speed, run_field_with_solver(timestep, GRAD_B, GET_SPEED, solve_diffeq_by_rkn), 1000, "4th order linearised Runge--Kutta--Nystrom (grad B, speed)");