Coryab/code #11

Merged
coryab merged 5 commits from coryab/code into develop 2023-10-24 18:35:52 +00:00
4 changed files with 163 additions and 341 deletions
Showing only changes of commit 6b6231fcf3 - Show all commits

View File

@ -46,7 +46,7 @@ instrument:
main: main.o $(LIBOBJS) $(CLASSOBJS) main: main.o $(LIBOBJS) $(CLASSOBJS)
$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
main_no_interface: main_no_interface.o $(LIBOBJS) $(CLASSOBJS) frequency_narrow_sweeps_long: frequency_narrow_sweeps_long.o $(LIBOBJS) $(CLASSOBJS)
$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)

View File

@ -0,0 +1,154 @@
/** @file main.cpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 1.0
*
* @brief The main program for this project
*
* @bug No known bugs
* */
#include <cmath>
#include <complex>
#include <fstream>
#include <omp.h>
#include <string>
#include <vector>
#include "PenningTrap.hpp"
#include "constants.hpp"
#include "utils.hpp"
#define PARTICLES 100
#define N 40000
/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time
* dependent potential.
*
* @details The simulation sweeps over different frequencies in [1., 1.7]
* MHz.
*
* */
void potential_resonance_narrow_sweep()
{
double time = 500.;
double amplitudes[]{.1, .4, .7};
double freq_start = 1.1;
double freq_end = 1.7;
double freq_increment = .0005;
size_t freq_iterations =
(size_t)((freq_end - freq_start) / freq_increment) + 1;
double res[4][freq_iterations];
std::string path = "output/time_dependent_potential/";
mkpath(path);
std::ofstream ofile;
#pragma omp parallel for
// Insert frequencies
for (size_t i = 0; i < freq_iterations; i++) {
res[0][i] = freq_start + freq_increment * i;
}
#pragma omp parallel
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((uint)PARTICLES);
#pragma omp for collapse(2)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {
// Reset particles and give new time dependent potential.
trap.reinitialize(amplitudes[i], res[0][j]);
res[i + 1][j] =
trap.fraction_of_particles_left(time, N, "rk4", false);
}
}
}
// Write results to file
ofile.open(path + "narrow_sweep_fine.txt");
for (size_t i = 0; i < freq_iterations; i++) {
ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ','
<< res[3][i] << '\n';
}
ofile.close();
}
/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time
* dependent potential.
*
* @details The simulation sweeps over different frequencies in [1., 1.7]
* MHz.
*
* */
void potential_resonance_narrow_sweep_interaction()
{
double time = 500.;
double amplitudes[]{.1, .4, .7};
double freq_start = 1.1;
double freq_end = 1.7;
double freq_increment = .0005;
size_t freq_iterations =
(size_t)((freq_end - freq_start) / freq_increment) + 1;
double res[4][freq_iterations];
std::string path = "output/time_dependent_potential/";
mkpath(path);
std::ofstream ofile;
#pragma omp parallel for
for (size_t i = 0; i < freq_iterations; i++) {
res[0][i] = freq_start + freq_increment * i;
}
#pragma omp parallel
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((uint)PARTICLES);
#pragma omp for collapse(2)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {
// Reset particles and give new time dependent potential.
trap.reinitialize(amplitudes[i], res[0][j]);
res[i + 1][j] = trap.fraction_of_particles_left(time, N);
}
}
}
// Write results to file
ofile.open(path + "narrow_sweep_interactions_fine.txt");
for (size_t i = 0; i < freq_iterations; i++) {
ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ','
<< res[3][i] << '\n';
}
ofile.close();
}
int main()
{
double start, end;
start = omp_get_wtime();
potential_resonance_narrow_sweep();
potential_resonance_narrow_sweep_interaction();
end = omp_get_wtime();
std::cout << "Time: " << end - start << " seconds" << std::endl;
return 0;
}

View File

@ -1,13 +1,15 @@
#!/bin/bash #!/bin/bash
#SBATCH --account=ec54 #SBATCH --account=ec54
#SBATCH --job-name=simple #SBATCH --job-name=particle-sim
#SBATCH --time=0-00:05:00 #SBATCH --time=0-01:30:00
#SBATCH --mem-per-cpu=4G #SBATCH --mem-per-cpu=2G
#SBATCH --cpus-per-task=8 #SBATCH --cpus-per-task=16
set -o errexit # Exit the script on any error set -o errexit # Exit the script on any error
set -o nounset # Treat any unset variables as an error set -o nounset # Treat any unset variables as an error
module --quiet purge # Reset the modules to the system default
srun ./main_no_interface module --quiet purge # Reset the modules to the system default
module load Armadillo/11.4.3-foss-2022b
srun ./frequency_narrow_sweeps_long

View File

@ -1,334 +0,0 @@
/** @file main.cpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 1.0
*
* @brief The main program for this project
*
* @bug No known bugs
* */
#include <cmath>
#include <complex>
#include <fstream>
#include <omp.h>
#include <string>
#include <vector>
#include "PenningTrap.hpp"
#include "constants.hpp"
#include "utils.hpp"
#define PARTICLES 100
#define N 40000
// Particles used for testing
Particle p1(vec3{20., 0., 20.}, vec3{0., 25., 0.}); ///< Particle 1
Particle p2(vec3{25., 25., 0.}, vec3{0., 40., 5.}); ///< Particle 2
/** @brief The analytical solution for particle p1
*
* @param t Time
*
* @return vec3
* */
vec3 analytical_solution_particle_1(double t)
{
double w_0 = T / CA_MASS;
double w_z2 = (50. * V / 1000.) / (CA_MASS * 500. * 500.);
double w_p = (w_0 + std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.;
double w_n = (w_0 - std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.;
double A_p = (25. + w_n * 20.) / (w_n - w_p);
double A_n = -(25. + w_p * 20.) / (w_n - w_p);
std::complex<double> f =
A_p * std::exp(std::complex<double>(0., -w_p * t))
+ A_n * std::exp(std::complex<double>(0., -w_n * t));
vec3 res{std::real(f), std::imag(f), 20. * std::cos(std::sqrt(w_z2) * t)};
return res;
}
/** @brief Simulate a single particle over the period of 50 \f$ \mu s \f$.
* */
void simulate_single_particle()
{
// Initialize trap with particle 1
PenningTrap trap(std::vector<Particle>{p1});
double time = 50.; // microseconds
// Simulate and write results to file
trap.write_simulation_to_dir("output/simulate_single_particle", time, N,
"rk4", false);
}
/** @brief Simulate 2 particles over the period of 50 \f$ \mu s \f$ with and
* without particle interactions.
* */
void simulate_two_particles()
{
// Initialize traps with particles
PenningTrap trap_no_interaction(std::vector<Particle>{p1, p2});
PenningTrap trap_with_interaction(std::vector<Particle>{p1, p2});
double time = 50.; // microseconds
// Simulate and write results to files
trap_no_interaction.write_simulation_to_dir(
"output/simulate_2_particles/no_interaction", time, N, "rk4", false);
trap_with_interaction.write_simulation_to_dir(
"output/simulate_2_particles/with_interaction", time, N);
}
/** @brief Simulate a single particle over 50 \f$ \mu s \f$ using different
* amount of steps and different methods.
* */
void simulate_single_particle_with_different_steps()
{
double time = 50.; // microseconds
std::ofstream ofile;
// Calculate relative error for RK4
std::string path = "output/relative_error/RK4/";
mkpath(path);
#pragma omp parallel for private(ofile)
for (int i = 0; i < 4; i++) {
int steps = 4000 * std::pow(2, i);
std::cout << steps << std::endl;
double dt = time / (double)steps;
ofile.open(path + std::to_string(steps) + "_steps.txt");
PenningTrap trap(std::vector<Particle>{p1});
simulation_t res = trap.simulate(time, steps, "rk4", false);
for (int i = 0; i < steps; i++) {
ofile << arma::norm(res.r_vecs[0][i]
- analytical_solution_particle_1(dt * i))
<< '\n';
}
ofile.close();
}
// Calculate relative error for forward Euler
path = "output/relative_error/euler/";
mkpath(path);
#pragma omp parallel for private(ofile)
for (int i = 0; i < 4; i++) {
int steps = 4000 * std::pow(2, i);
double dt = time / (double)steps;
ofile.open(path + std::to_string(steps) + "_steps.txt");
PenningTrap trap(std::vector<Particle>{p1});
simulation_t res = trap.simulate(time, steps, "euler", false);
for (int i = 0; i < steps; i++) {
ofile << arma::norm(res.r_vecs[0][i]
- analytical_solution_particle_1(dt * i))
<< '\n';
}
ofile.close();
}
}
/** @brief Simulate 100 particles over 50 \f$ \mu s \f$.
* */
void simulate_100_particles()
{
PenningTrap trap((unsigned)100);
double time = 50.; // microseconds
// trap.write_simulation_to_dir("output/simulate_100_particles", time, N,
//"rk4", false);
trap.simulate(time, N, "rk4", true);
}
/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time
* dependent potential.
*
* @details The simulation sweeps over different frequencies in [0.2, 2.5]
* MHz.
*
* */
void potential_resonance_wide_sweep()
{
double time = 500.;
double amplitudes[]{.1, .4, .7};
double freq_start = .2;
double freq_end = 2.5;
double freq_increment = .02;
size_t freq_iterations =
(size_t)((freq_end - freq_start) / freq_increment) + 1;
double res[4][freq_iterations];
std::string path = "output/time_dependent_potential/";
mkpath(path);
std::ofstream ofile;
#pragma omp parallel for
// Insert frequencies
for (size_t i = 0; i < freq_iterations; i++) {
res[0][i] = freq_start + freq_increment * i;
}
#pragma omp parallel
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((uint)100);
#pragma omp for collapse(2)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {
// Reset particles and give new time dependent potential.
trap.reinitialize(amplitudes[i], res[0][j]);
res[i + 1][j] =
trap.fraction_of_particles_left(time, N, "rk4", false);
}
}
}
// Write results to file
ofile.open(path + "wide_sweep.txt");
for (size_t i = 0; i < freq_iterations; i++) {
ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ','
<< res[3][i] << '\n';
}
ofile.close();
}
/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time
* dependent potential.
*
* @details The simulation sweeps over different frequencies in [1., 1.7]
* MHz.
*
* */
void potential_resonance_narrow_sweep()
{
double time = 500.;
double amplitudes[]{.1, .4, .7};
double freq_start = 1.;
double freq_end = 1.7;
double freq_increment = .002;
size_t freq_iterations =
(size_t)((freq_end - freq_start) / freq_increment) + 1;
double res[4][freq_iterations];
std::string path = "output/time_dependent_potential/";
mkpath(path);
std::ofstream ofile;
#pragma omp parallel for
// Insert frequencies
for (size_t i = 0; i < freq_iterations; i++) {
res[0][i] = freq_start + freq_increment * i;
}
#pragma omp parallel
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((uint)100);
#pragma omp for collapse(2)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {
// Reset particles and give new time dependent potential.
trap.reinitialize(amplitudes[i], res[0][j]);
res[i + 1][j] =
trap.fraction_of_particles_left(time, N, "rk4", false);
}
}
}
// Write results to file
ofile.open(path + "narrow_sweep.txt");
for (size_t i = 0; i < freq_iterations; i++) {
ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ','
<< res[3][i] << '\n';
}
ofile.close();
}
/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time
* dependent potential.
*
* @details The simulation sweeps over different frequencies in [1., 1.7]
* MHz.
*
* */
void potential_resonance_narrow_sweep_interaction()
{
double time = 500.;
double amplitudes[]{.1, .4, .7};
double freq_start = 1.;
double freq_end = 1.7;
double freq_increment = .002;
size_t freq_iterations =
(size_t)((freq_end - freq_start) / freq_increment) + 1;
double res[4][freq_iterations];
std::string path = "output/time_dependent_potential/";
mkpath(path);
std::ofstream ofile;
#pragma omp parallel for
for (size_t i = 0; i < freq_iterations; i++) {
res[0][i] = freq_start + freq_increment * i;
}
#pragma omp parallel
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((uint)100);
#pragma omp for collapse(2)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {
// Reset particles and give new time dependent potential.
trap.reinitialize(amplitudes[i], res[0][j]);
res[i + 1][j] = trap.fraction_of_particles_left(time, N);
}
}
}
// Write results to file
ofile.open(path + "narrow_sweep_interactions.txt");
for (size_t i = 0; i < freq_iterations; i++) {
ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ','
<< res[3][i] << '\n';
}
ofile.close();
}
int main()
{
int option = 1;
bool chosen = false;
double start, end;
start = omp_get_wtime();
simulate_100_particles();
potential_resonance_wide_sweep();
// potential_resonance_narrow_sweep();
// potential_resonance_narrow_sweep_interaction();
end = omp_get_wtime();
std::cout << "Time: " << end - start << " seconds" << std::endl;
return 0;
}