Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • engelerp/rbcomb-simulation
1 result
Show changes
Commits on Source (8)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1156 additions and 855 deletions
parameter_calculation/Untitled.ipynb
View file @ 5f273894
...@@ -328,7 +328,7 @@ ...@@ -328,7 +328,7 @@
"name": "python", "name": "python",
"nbconvert_exporter": "python", "nbconvert_exporter": "python",
"pygments_lexer": "ipython3", "pygments_lexer": "ipython3",
"version": "3.7.4" "version": "3.9.4"
} }
}, },
"nbformat": 4, "nbformat": 4,
......
projects/braidingTightBinding/bin/Makefile
View file @ 5f273894
...@@ -2,7 +2,7 @@ INCLUDEDIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/inclu ...@@ -2,7 +2,7 @@ INCLUDEDIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/inclu
LIBDIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/lib LIBDIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/lib
HDF5DIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/lib/src HDF5DIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/lib/src
CXX = /usr/local/bin/g++-9 CXX = /usr/local/bin/g++-10
CXXFLAGS = -std=c++2a -Wall -Wpedantic -I$(INCLUDEDIR) -I$(LIBDIR) -I${HDF5DIR} -O5 -mtune=native -march=native -framework Accelerate -lblas CXXFLAGS = -std=c++2a -Wall -Wpedantic -I$(INCLUDEDIR) -I$(LIBDIR) -I${HDF5DIR} -O5 -mtune=native -march=native -framework Accelerate -lblas
run: all run: all
...@@ -11,10 +11,22 @@ run: all ...@@ -11,10 +11,22 @@ run: all
./main_braiding_zerostate ../results/data/zerostate/ ./main_braiding_zerostate ../results/data/zerostate/
python3 plot_zerostate.py python3 plot_zerostate.py
all: main_braiding_zerostate run_nocoupling: main_nocoupling
mkdir -p ../results/data/nocoupling/
mkdir -p ../results/plots/nocoupling/
./main_nocoupling ../results/data/nocoupling/
python3 plot.py ../results/data/nocoupling/ ../results/plots/nocoupling/ "No Coupling" 0 10
main_braiding_zerostate: main_braiding_zerostate.cpp ${LIBDIR}/diagonalizer.hpp ${LIBDIR}/drum.hpp ${INCLUDEDIR}/grabber.hpp ${LIBDIR}/rk4_buffer.hpp ${LIBDIR}/rk4_stepper.hpp ${LIBDIR}/system.hpp ${LIBDIR}/system_parameters.hpp ${LIBDIR}/vec2.hpp ${INCLUDEDIR}/coupler_braid.hpp ${INCLUDEDIR}/driver_braid.hpp ${INCLUDEDIR}/drum_parameters_braid.hpp ${INCLUDEDIR}/drum_variables_braid.hpp ${INCLUDEDIR}/force_braid.hpp ${INCLUDEDIR}/matrix_element_calculator_braid.hpp ${INCLUDEDIR}/rbcomb_generator_braid.hpp ${INCLUDEDIR}/vortex.hpp all: main_braiding_zerostate main_nocoupling
main_braiding_zerostate: main_braiding_zerostate.cpp ${LIBDIR}/diagonalizer.hpp ${LIBDIR}/drum.hpp ${INCLUDEDIR}/grabber.hpp ${LIBDIR}/rk4_buffer.hpp ${LIBDIR}/rk4_stepper.hpp ${LIBDIR}/system.hpp ${LIBDIR}/system_parameters.hpp ${LIBDIR}/vec2.hpp ${INCLUDEDIR}/coupler_braid.hpp ${INCLUDEDIR}/coupler_const.hpp ${INCLUDEDIR}/driver_braid.hpp ${INCLUDEDIR}/drum_parameters_braid.hpp ${INCLUDEDIR}/drum_variables_braid.hpp ${INCLUDEDIR}/force_braid.hpp ${INCLUDEDIR}/matrix_element_calculator_braid.hpp ${INCLUDEDIR}/rbcomb_generator_braid.hpp ${INCLUDEDIR}/vortex.hpp
$(CXX) $< $(CXXFLAGS) -o $@
main_nocoupling: main_nocoupling.cpp ${LIBDIR}/diagonalizer.hpp ${LIBDIR}/drum.hpp ${INCLUDEDIR}/grabber.hpp ${LIBDIR}/rk4_buffer.hpp ${LIBDIR}/rk4_stepper.hpp ${LIBDIR}/system.hpp ${LIBDIR}/system_parameters.hpp ${LIBDIR}/vec2.hpp ${INCLUDEDIR}/coupler_const.hpp ${INCLUDEDIR}/driver_braid.hpp ${INCLUDEDIR}/drum_parameters_braid.hpp ${INCLUDEDIR}/drum_variables_braid.hpp ${INCLUDEDIR}/force_braid.hpp ${INCLUDEDIR}/matrix_element_calculator_braid.hpp ${INCLUDEDIR}/rbcomb_generator_braid.hpp
$(CXX) $< $(CXXFLAGS) -o $@ $(CXX) $< $(CXXFLAGS) -o $@
plot:
python3 plot_zerostate.py
clean: clean:
rm -rf main_braiding_zerostate ../results/data/* ../results/plots/* rm -rf main_braiding_zerostate ../results/data/* ../results/plots/*
projects/braidingTightBinding/bin/main_braiding_zerostate deleted 100755 → 0
View file @ 92fc03a4
File deleted
projects/braidingTightBinding/bin/main_braiding_zerostate.cpp
View file @ 5f273894
...@@ -25,111 +25,121 @@ ...@@ -25,111 +25,121 @@
template <typename iterator_t> template <typename iterator_t>
void writeToFile(const std::string filename, const std::string header, int line_length, iterator_t begin, iterator_t end) noexcept void writeToFile(const std::string filename, const std::string header, int line_length, iterator_t begin, iterator_t end) noexcept
{ {
//no line breaks if line_length == -1 //no line breaks if line_length == -1
if(line_length == -1) if(line_length == -1){
line_length = std::distance(begin, end); line_length = std::distance(begin, end);
std::fstream file (filename, file.out); }
file << header; std::fstream file (filename, file.out);
file << std::endl; file << header;
for(auto it = begin; it != end; ++it){ file << std::endl;
if(std::distance(begin, it) % line_length == 0) for(auto it = begin; it != end; ++it){
file << std::endl; if(std::distance(begin, it) % line_length == 0){
file << *it << " "; file << std::endl;
} }
file << std::flush; file << *it << " ";
file.close(); }
file << std::flush;
file.close();
} }
int main(int argc, char** argv){ int main(int argc, char** argv){
assert(argc == 2); assert(argc == 2);
//shorthands //shorthands
using value_t = double; using value_t = double;
using params_t = DrumParametersBraid<value_t>; using params_t = DrumParametersBraid<value_t>;
using vars_t = DrumVariablesBraid<value_t>; using vars_t = DrumVariablesBraid<value_t>;
using buffer_t = RK4Buffer<value_t>; using buffer_t = RK4Buffer<value_t>;
using drum_t = Drum<value_t, params_t, vars_t, buffer_t>; using drum_t = Drum<value_t, params_t, vars_t, buffer_t>;
using coupler_t = CouplerBraid<value_t, drum_t, params_t>; using coupler_t = CouplerBraid<value_t, drum_t, params_t>;
using driver_t = DriverBraid<value_t, drum_t>; using driver_t = DriverBraid<value_t, drum_t>;
using sysparams_t = SystemParameters<coupler_t, driver_t>; using sysparams_t = SystemParameters<coupler_t, driver_t>;
using force_t = ForceBraid<value_t, params_t, vars_t, buffer_t>; using force_t = ForceBraid<value_t, params_t, vars_t, buffer_t>;
using stepper_t = Rk4Stepper<value_t, params_t, vars_t, buffer_t, force_t>; using stepper_t = Rk4Stepper<value_t, params_t, vars_t, buffer_t, force_t>;
using generator_t = RbcombGeneratorBraid<value_t, params_t, vars_t, buffer_t>; using generator_t = RbcombGeneratorBraid<value_t, params_t, vars_t, buffer_t>;
using grabber_t = Grabber<value_t, drum_t>; using grabber_t = Grabber<value_t, drum_t>;
using matelecalc_t = MatrixElementCalculatorBraid<value_t, params_t, vars_t, drum_t>; using matelecalc_t = MatrixElementCalculatorBraid<value_t, params_t, vars_t, drum_t>;
using system_t = System<value_t, drum_t, grabber_t, sysparams_t, force_t, coupler_t, driver_t, stepper_t, matelecalc_t>; using system_t = System<value_t, drum_t, grabber_t, sysparams_t, force_t, coupler_t, driver_t, stepper_t, matelecalc_t>;
using vortex_t = Vortex<value_t>; using vortex_t = Vortex<value_t>;
using vec2_t = Vec2<value_t>; using vec2_t = Vec2<value_t>;
value_t k0 =24674011002.723392; value_t k0 = 24674011002.723392;
value_t k1 = 493480220.0544679; value_t k1 = 493480220.0544679;
value_t k2 = 2467401100.2723393; value_t k2 = 2467401100.2723393;
value_t c = 5.235987755982988; value_t c = 5.235987755982988;
vortex_t vortex0 (vec2_t(10.5,24.), 5., 1., 1., 1.); //as Eliska did it vortex_t vortex0 (vec2_t(10.5,24.), 5., 1., 1., 1.); //as Eliska did it
vortex_t vortex1 (vec2_t(31.5,24.), 5., 1., 1., 1.); //as Eliska did it vortex_t vortex1 (vec2_t(31.5,24.), 5., 1., 1., 1.); //as Eliska did it
vortex_t vortex2 (vec2_t(52.5,24.), 5., 1., 1., 1.); //as Eliska did it vortex_t vortex2 (vec2_t(52.5,24.), 5., 1., 1., 1.); //as Eliska did it
//std::vector<vortex_t> vortices ({vortex0, vortex1, vortex2}); //std::vector<vortex_t> vortices ({vortex0, vortex1, vortex2});
std::vector<vortex_t> vortices ({vortex0, vortex1, vortex2}); std::vector<vortex_t> vortices ({vortex1});
params_t drum_parameters (k0, k1, k2, c, vec2_t (0., 0.), 'A'); params_t drum_parameters (k0, k1, k2, c, vec2_t (0., 0.), 'A');
generator_t lattice_generator (37,15); generator_t lattice_generator (37,15);
auto lat_pair = lattice_generator(drum_parameters); auto lat_pair = lattice_generator(drum_parameters);
//change edge potentials to push edge modes around //change edge potentials to push edge modes around
value_t onsite_A = 27388152213.022964; value_t onsite_A = 27388152213.022964;
value_t onsite_B = 21959869792.42382; value_t onsite_B = 21959869792.42382;
for(size_t i = 0; i < lat_pair.first.size()-1; ++i){ for(size_t i = 0; i < lat_pair.first.size()-1; ++i){
if(lat_pair.second[i][0] == lat_pair.first.size() || lat_pair.second[i][1] == lat_pair.first.size() || lat_pair.second[i][2] == lat_pair.first.size()){ if(lat_pair.second[i][0] == lat_pair.first.size() || lat_pair.second[i][1] == lat_pair.first.size() || lat_pair.second[i][2] == lat_pair.first.size()){
if(lat_pair.first[i].get_parameters().sublattice == 'A'){ //'A' site if(lat_pair.first[i].get_parameters().sublattice == 'A'){ //'A' site
lat_pair.first[i].get_parameters().k0 = onsite_A; lat_pair.first[i].get_parameters().k0 = onsite_A;
for(int j: lat_pair.second[i]) for(int j: lat_pair.second[i]){
if(j != lat_pair.first.size()) if(j != lat_pair.first.size()){
lat_pair.first[j].get_parameters().k0 = onsite_B; lat_pair.first[j].get_parameters().k0 = onsite_B;
}
}
} }
else{ //'B' site else{ //'B' site
lat_pair.first[i].get_parameters().k0 = onsite_B; lat_pair.first[i].get_parameters().k0 = onsite_B;
for(int j: lat_pair.second[i]) for(int j: lat_pair.second[i]){
if(j != lat_pair.first.size()) if(j != lat_pair.first.size()){
lat_pair.first[j].get_parameters().k0 = onsite_A; lat_pair.first[j].get_parameters().k0 = onsite_A;
}
}
} }
} }
} }
//set initial conditions
double x_0 = 0.;
double v_0 = 0.;
for(auto drum: lat_pair.first){
drum.get_variables().x = x_0;
drum.get_variables().xdot = v_0;
}
sysparams_t system_parameters (coupler_t(vortices), driver_t(), lat_pair.second);
force_t force;
stepper_t stepper;
//generate grabber sysparams_t system_parameters (coupler_t(vortices), driver_t(), lat_pair.second);
size_t save_every = 10; force_t force;
std::string parameters_file = std::string(argv[1]) + "PARAMS.TXT"; stepper_t stepper;
std::string adjacency_file = std::string(argv[1]) + "ADJACENCY.TXT";
std::string dynamics_file = std::string(argv[1]) + "DYNAMICS.TXT";
grabber_t grabber (save_every, parameters_file, adjacency_file, dynamics_file);
system_t system (1., 0.00001, lat_pair.first, stepper, force, system_parameters, grabber); //generate grabber
system.save(); size_t save_every = 10;
std::string parameters_file = std::string(argv[1]) + "PARAMS.TXT";
std::string adjacency_file = std::string(argv[1]) + "ADJACENCY.TXT";
std::string dynamics_file = std::string(argv[1]) + "DYNAMICS.TXT";
grabber_t grabber (save_every, parameters_file, adjacency_file, dynamics_file);
//set up matrix element calculator system_t system (1., 0.00001, lat_pair.first, stepper, force, system_parameters, grabber);
matelecalc_t mec; system.save();
//set up diagonalizer //set up matrix element calculator
Diagonalizer diagonalizer; matelecalc_t mec;
//calculate eigenvalues
std::vector<value_t> matrix = system.get_matrix(mec);
//std::vector<value_t> evals = diagonalizer.ev(system.get_matrix(mec), lat_pair.first.size()-1);
std::pair<std::vector<value_t>, std::vector<value_t> > diag_pair = diagonalizer.evv(system.get_matrix(mec), lat_pair.first.size()-1);
//make eigenvalues to frequencies //set up diagonalizer
std::transform(diag_pair.first.begin(), diag_pair.first.end(), diag_pair.first.begin(), [](value_t x) -> value_t {return std::sqrt(x)/(2.*M_PI);}); Diagonalizer diagonalizer;
//calculate eigenvalues
std::vector<value_t> matrix = system.get_matrix(mec);
//std::vector<value_t> evals = diagonalizer.ev(system.get_matrix(mec), lat_pair.first.size()-1);
std::pair<std::vector<value_t>, std::vector<value_t> > diag_pair = diagonalizer.evv(system.get_matrix(mec), lat_pair.first.size()-1);
//find the interesting eigenvector range //make eigenvalues to frequencies
size_t ev_begin = 913 * diag_pair.first.size(); std::transform(diag_pair.first.begin(), diag_pair.first.end(), diag_pair.first.begin(), [](value_t x) -> value_t {return std::sqrt(x)/(2.*M_PI);});
size_t ev_end = 943 * diag_pair.first.size();
//write to file //write to file
writeToFile(std::string(argv[1]) + "eigenvalues.txt", "", 1, diag_pair.first.begin(), diag_pair.first.end()); writeToFile(std::string(argv[1]) + "eigenvalues.txt", "", 1, diag_pair.first.begin(), diag_pair.first.end());
writeToFile(std::string(argv[1]) + "eigenvectors.txt", "", diag_pair.first.size(), diag_pair.second.begin(), diag_pair.second.end()); writeToFile(std::string(argv[1]) + "eigenvectors.txt", "", diag_pair.first.size(), diag_pair.second.begin(), diag_pair.second.end());
} }
projects/braidingTightBinding/bin/main_nocoupling.cpp 0 → 100644
View file @ 5f273894
#include <coupler.hpp>
#include <diagonalizer.hpp>
#include <driver.hpp>
#include <drum.hpp>
#include <force.hpp>
#include <grabber.hpp>
#include <lattice_generator.hpp>
#include <rk4_buffer.hpp>
#include <rk4_stepper.hpp>
#include <system.hpp>
#include <system_parameters.hpp>
#include <vec2.hpp>
#include <coupler_const.hpp>
#include <driver_braid.hpp>
#include <drum_parameters_braid.hpp>
#include <drum_variables_braid.hpp>
#include <force_braid.hpp>
#include <matrix_element_calculator_braid.hpp>
#include <rbcomb_generator_braid.hpp>
#include <vortex.hpp>
#include <string>
#include <fstream>
#include <cassert>
template <typename iterator_t>
void writeToFile(const std::string filename, const std::string header, int line_length, iterator_t begin, iterator_t end) noexcept
{
//no line breaks if line_length == -1
if(line_length == -1){
line_length = std::distance(begin, end);
}
std::fstream file (filename, file.out);
file << header;
file << std::endl;
for(auto it = begin; it != end; ++it){
if(std::distance(begin, it) % line_length == 0){
file << std::endl;
}
file << *it << " ";
}
file << std::flush;
file.close();
}
int main(int argc, char** argv){
assert(argc == 2);
//shorthands
using value_t = double;
using params_t = DrumParametersBraid<value_t>;
using vars_t = DrumVariablesBraid<value_t>;
using buffer_t = RK4Buffer<value_t>;
using drum_t = Drum<value_t, params_t, vars_t, buffer_t>;
using coupler_t = CouplerConst<value_t, drum_t, params_t>;
using driver_t = DriverBraid<value_t, drum_t>;
using sysparams_t = SystemParameters<coupler_t, driver_t>;
using force_t = ForceBraid<value_t, params_t, vars_t, buffer_t>;
using stepper_t = Rk4Stepper<value_t, params_t, vars_t, buffer_t, force_t>;
using generator_t = RbcombGeneratorBraid<value_t, params_t, vars_t, buffer_t>;
using grabber_t = Grabber<value_t, drum_t>;
using matelecalc_t = MatrixElementCalculatorBraid<value_t, params_t, vars_t, drum_t>;
using system_t = System<value_t, drum_t, grabber_t, sysparams_t, force_t, coupler_t, driver_t, stepper_t, matelecalc_t>;
using vortex_t = Vortex<value_t>;
using vec2_t = Vec2<value_t>;
value_t k0 = 24674011002.723392;
value_t k1 = 0.;
value_t k2 = 2467401100.2723393;
value_t c = 5.235987755982988;
params_t drum_parameters (k0, k1, k2, c, vec2_t (0., 0.), 'A');
generator_t lattice_generator (37,15);
auto lat_pair = lattice_generator(drum_parameters);
//set initial conditions
double x_0 = 1.;
double xdot_0 = 0.;
for(drum_t& drum: lat_pair.first){
drum.get_variables().x = x_0;
drum.get_variables().xdot = xdot_0;
}
sysparams_t system_parameters (coupler_t(0.), driver_t(), lat_pair.second);
force_t force;
stepper_t stepper;
//generate grabber
double period = 0.00004;
double points_per_period = 100.;
size_t save_every = 10;
std::string parameters_file = std::string(argv[1]) + "PARAMS.TXT";
std::string adjacency_file = std::string(argv[1]) + "ADJACENCY.TXT";
std::string dynamics_file = std::string(argv[1]) + "DYNAMICS.TXT";
grabber_t grabber (save_every, parameters_file, adjacency_file, dynamics_file);
system_t system (10000*period, period/points_per_period, lat_pair.first, stepper, force, system_parameters, grabber);
system.simulate();
system.save();
//set up matrix element calculator
matelecalc_t mec;
//set up diagonalizer
Diagonalizer diagonalizer;
//calculate eigenvalues
std::vector<value_t> matrix = system.get_matrix(mec);
//std::vector<value_t> evals = diagonalizer.ev(system.get_matrix(mec), lat_pair.first.size()-1);
std::pair<std::vector<value_t>, std::vector<value_t> > diag_pair = diagonalizer.evv(system.get_matrix(mec), lat_pair.first.size()-1);
//make eigenvalues to frequencies
std::transform(diag_pair.first.begin(), diag_pair.first.end(), diag_pair.first.begin(), [](value_t x) -> value_t {return std::sqrt(x)/(2.*M_PI);});
//write to file
writeToFile(std::string(argv[1]) + "eigenvalues.txt", "", 1, diag_pair.first.begin(), diag_pair.first.end());
writeToFile(std::string(argv[1]) + "eigenvectors.txt", "", diag_pair.first.size(), diag_pair.second.begin(), diag_pair.second.end());
}
projects/braidingTightBinding/bin/plot.py 0 → 100644
View file @ 5f273894
import matplotlib.pyplot as plt
import numpy as np
import sys
if len(sys.argv) < 4:
print(f"Usage: python3.9 {argv[0]} ../path/to/data/ ../path/to/plots/ name [modes_range_start modes_range_end]")
sys.exit()
readpath = sys.argv[1]
writepath = sys.argv[2]
name = sys.argv[3]
#eigenvalues
ev = np.loadtxt(readpath+"eigenvalues.txt")
starting_mode = 1110
end_mode = 1180
if len(sys.argv) == 6:
starting_mode = int(sys.argv[4])
end_mode = int(sys.argv[5])
if end_mode - starting_mode > 70:
print(f"Maximum number of modes: 70 (entered: {end_mode - starting_mode})")
end_mode = starting_mode + 70
print(f"Using Modes {starting_mode}-{end_mode}")
plt.plot(np.array(list(range(starting_mode,end_mode))),ev[starting_mode:end_mode],'o', ms=1)
plt.title(f"Frequency levels under {name}")
plt.xlabel("Level index")
plt.ylabel("Frequency [Hz]")
plt.savefig(writepath + "Frequencies_zoomed.pdf", dpi = 500)
plt.clf()
plt.plot(ev[:],'o', ms=1)
plt.title(f"Frequency levels under {name}")
plt.xlabel("Level index")
plt.ylabel("Frequency [Hz]")
plt.savefig(writepath + "Frequencies.pdf", dpi = 500)
#eigenvectors
data = np.loadtxt(readpath+"eigenvectors.txt")
pos_data = np.loadtxt(readpath+"PARAMS.txt", dtype=str) #positions are rows 4 and 5
desired_length = len(data[0,:])
pos_x = pos_data[0:desired_length,4].astype('float')
pos_y = pos_data[:desired_length,5].astype('float')
plt.figure(figsize=(8,5*(end_mode - starting_mode + 1)))
for i in range(starting_mode,end_mode):
plt.subplot(int(end_mode-starting_mode),1,i-starting_mode+1)
cmax = np.amax(data[i,:])
cmin = np.amin(data[i,:])
if(-cmin > cmax):
cmax = -cmin
elif(-cmin <= cmax):
cmin = -cmax
plt.scatter(pos_x, pos_y, c=data[i,:], edgecolors='k', linewidth=0.1, cmap='seismic', s=16, vmin=cmin, vmax=cmax)
plt.colorbar()
plt.title(f"Mode {i+1}, ({ev[i]} Hz)")
plt.xlabel("x")
plt.ylabel("y")
plt.savefig(writepath + f"Modes.pdf", dpi = 40, bbox_inches='tight', pad_inches=0)
projects/braidingTightBinding/bin/plot_states.py 0 → 100644
View file @ 5f273894
import matplotlib.pyplot as plt
import numpy as np
import sys
readpath = "../results/data/zerostate/"
writepath = "../results/plots/zerostate/"
if(len(sys.argv) < 2):
print(f"Usage: python3.9 {sys.argv[0]} mode1 mode2 mode3 ...")
else:
num_modes = int(len(sys.argv)-1)
#eigenvalues
ev = np.loadtxt(readpath+"eigenvalues.txt")
#eigenvectors
data = np.loadtxt(readpath+"eigenvectors.txt")
pos_data = np.loadtxt(readpath+"PARAMS.txt", dtype=str) #positions are rows 4 and 5
desired_length = len(data[0,:])
pos_x = pos_data[0:desired_length,4].astype('float')
pos_y = pos_data[:desired_length,5].astype('float')
plt.figure(figsize=(8,5*num_modes))
count = 1
for str_i in sys.argv[1:]:
i = int(str_i)-1
plt.subplot(num_modes, 1, count)
cmax = np.amax(data[i,:])
cmin = np.amin(data[i,:])
if(-cmin > cmax):
cmax = -cmin
elif(-cmin <= cmax):
cmin = -cmax
plt.scatter(pos_x, pos_y, c=data[i,:], edgecolors='k', linewidth=0.1, cmap='seismic', s=16, vmin=cmin, vmax=cmax)
plt.colorbar()
plt.title(f"Mode {i+1}, ({ev[i]} Hz)")
plt.xlabel("x")
plt.ylabel("y")
count += 1
name = "Modes"
for str_i in sys.argv[1:]:
name += f"_{str_i}"
plt.savefig(writepath + name, dpi=500, bbox_inches='tight', pad_inches=0)
projects/braidingTightBinding/bin/plot_zerostate.py
View file @ 5f273894
...@@ -7,8 +7,8 @@ writepath = "../results/plots/zerostate/" ...@@ -7,8 +7,8 @@ writepath = "../results/plots/zerostate/"
#eigenvalues #eigenvalues
ev = np.loadtxt(readpath+"eigenvalues.txt") ev = np.loadtxt(readpath+"eigenvalues.txt")
starting_mode = 1115 starting_mode = 1110
end_mode = 1174 end_mode = 1180
plt.plot(np.array(list(range(starting_mode,end_mode))),ev[starting_mode:end_mode],'o', ms=1) plt.plot(np.array(list(range(starting_mode,end_mode))),ev[starting_mode:end_mode],'o', ms=1)
plt.title("Frequency levels under vortex Kekule") plt.title("Frequency levels under vortex Kekule")
...@@ -35,7 +35,7 @@ pos_y = pos_data[:desired_length,5].astype('float') ...@@ -35,7 +35,7 @@ pos_y = pos_data[:desired_length,5].astype('float')
plt.figure(figsize=(8,300)) plt.figure(figsize=(8,300))
for i in range(starting_mode,end_mode): for i in range(starting_mode,end_mode):
plt.subplot(60,1,i-starting_mode+1) plt.subplot(int(end_mode-starting_mode),1,i-starting_mode+1)
cmax = np.amax(data[i,:]) cmax = np.amax(data[i,:])
cmin = np.amin(data[i,:]) cmin = np.amin(data[i,:])
if(-cmin > cmax): if(-cmin > cmax):
...@@ -48,4 +48,4 @@ for i in range(starting_mode,end_mode): ...@@ -48,4 +48,4 @@ for i in range(starting_mode,end_mode):
plt.xlabel("x") plt.xlabel("x")
plt.ylabel("y") plt.ylabel("y")
plt.savefig(writepath + f"Modes.pdf", dpi = 50, bbox_inches='tight', pad_inches=0) plt.savefig(writepath + f"Modes.pdf", dpi = 40, bbox_inches='tight', pad_inches=0)
projects/braidingTightBinding/include/coupler_braid.hpp
View file @ 5f273894
...@@ -6,55 +6,57 @@ ...@@ -6,55 +6,57 @@
template <typename value_t, typename drum_t, typename params_t> template <typename value_t, typename drum_t, typename params_t>
class CouplerBraid: public Coupler<value_t, drum_t>{ class CouplerBraid: public Coupler<value_t, drum_t>{
public: public:
//constructors //constructors
CouplerBraid(const std::vector<Vortex<value_t> >& vortices): vortices_(vortices) {} CouplerBraid(const std::vector<Vortex<value_t> >& vortices): vortices_(vortices) {}
~CouplerBraid() = default; ~CouplerBraid() = default;
void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
{ {
//we need the drum parameters later on, so extract it //we need the drum parameters later on, so extract it
params_.reserve(drum_vec.size()); params_.reserve(drum_vec.size());
for(auto d: drum_vec) for(auto d: drum_vec){
params_.push_back(d.get_parameters()); params_.push_back(d.get_parameters());
}
//we also need the neighbour vectors nicely indexed
for(auto d: drum_vec){ //we also need the neighbour vectors nicely indexed
if(d.get_parameters().sublattice == 'A'){ for(auto d: drum_vec){
n_vecs_.push_back(d.get_parameters().s0); if(d.get_parameters().sublattice == 'A'){
n_vecs_.push_back(d.get_parameters().s1); n_vecs_.push_back(d.get_parameters().s0);
n_vecs_.push_back(d.get_parameters().s2); n_vecs_.push_back(d.get_parameters().s1);
break; n_vecs_.push_back(d.get_parameters().s2);
} break;
} }
} }
void step(value_t dt) noexcept final }
{ void step(value_t dt) noexcept final
//move vortices or alpha here {
} //move vortices or alpha here
value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override }
{ value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
Vec2<value_t> v = params_[drum_index].position; {
Vec2<value_t> s = n_vecs_[neighbour_index]; Vec2<value_t> v = params_[drum_index].position;
Vec2<value_t> s = n_vecs_[neighbour_index];
//the formula is only valid on the 'A' sublattice
if(params_[drum_index].sublattice == 'B') //the formula is only valid on the 'A' sublattice
v = v - s; if(params_[drum_index].sublattice == 'B'){
v = v - s;
std::complex<value_t> ret_val = 1.; //vortices are multiplicative }
for(auto tex: vortices_){
ret_val *= tex.distortion(v); std::complex<value_t> ret_val = 1.; //vortices are multiplicative
} for(auto tex: vortices_){
//add kekule ret_val *= tex.distortion(v);
ret_val *= vortices_.begin()->kekule(v, s); }
//add kekule
return 1. + std::real<value_t>(ret_val); //return with added offset, as in eliska's code ret_val *= vortices_.begin()->kekule(v, s);
}
return 1. + std::real<value_t>(ret_val); //return with added offset, as in eliska's code
private: }
std::vector<Vortex<value_t> > vortices_;
std::vector<params_t> params_; private:
std::vector<Vec2<value_t> > n_vecs_; std::vector<Vortex<value_t> > vortices_;
std::vector<params_t> params_;
std::vector<Vec2<value_t> > n_vecs_;
}; };
#endif #endif
projects/braidingTightBinding/include/coupler_const.hpp 0 → 100644
View file @ 5f273894
#ifndef COUPLER_CONST_HPP_INCLUDED
#define COUPLER_CONST_HPP_INCLUDED
#include <vector>
#include <coupler.hpp>
template <typename value_t, typename drum_t, typename params_t>
class CouplerConst: public Coupler<value_t, drum_t>{
public:
//constructors
CouplerConst(value_t value): value_(value) {}
~CouplerConst() = default;
void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
{
//nothing to do
}
void step(value_t dt) noexcept final
{
//nothing to do
}
value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
{
return value_;
}
private:
value_t value_;
};
#endif
projects/braidingTightBinding/include/driver_braid.hpp
View file @ 5f273894
...@@ -5,23 +5,23 @@ ...@@ -5,23 +5,23 @@
template <typename value_t, typename drum_t> template <typename value_t, typename drum_t>
class DriverBraid: public Driver<value_t, drum_t>{ class DriverBraid: public Driver<value_t, drum_t>{
public: public:
//constructors //constructors
DriverBraid() = default; DriverBraid() = default;
~DriverBraid() = default; ~DriverBraid() = default;
void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
{ {
} }
void step(value_t dt) noexcept final void step(value_t dt) noexcept final
{ {
} }
value_t operator()(const size_t drum_index) const noexcept final override value_t operator()(const size_t drum_index) const noexcept final override
{ {
return 0; return 0;
} }
}; };
#endif #endif
projects/braidingTightBinding/include/drum_parameters_braid.hpp
View file @ 5f273894
...@@ -7,58 +7,59 @@ ...@@ -7,58 +7,59 @@
//TODO: Create an interface. It would be a trivial one, though. //TODO: Create an interface. It would be a trivial one, though.
template <typename value_t> template <typename value_t>
class DrumParametersBraid{ class DrumParametersBraid{
public: public:
/* Arguments: /* Arguments:
* k0: diagonal matrix element * k0: diagonal matrix element
* k1: coupling matrix element prefactor 1 * k1: coupling matrix element prefactor 1
* k2: coupling matrix element prefactor 2 * k2: coupling matrix element prefactor 2
* c: damping coefficient * c: damping coefficient
* position: position of the drum in real space * position: position of the drum in real space
* sublattice: 'A'/'a' or 'B'/'b', identifies sublattice drum is part of * sublattice: 'A'/'a' or 'B'/'b', identifies sublattice drum is part of
*/ */
DrumParametersBraid(const value_t k0, const value_t k1, const value_t k2, const value_t c, const Vec2<value_t>& position, char sublattice) noexcept DrumParametersBraid(const value_t k0, const value_t k1, const value_t k2, const value_t c, const Vec2<value_t>& position, char sublattice) noexcept
: k0(k0), k1(k1), k2(k2), c(c), position(position), sublattice(sublattice) : k0(k0), k1(k1), k2(k2), c(c), position(position), sublattice(sublattice)
{ {
generate_sublattice(sublattice, value_t(1.)); generate_sublattice(sublattice, value_t(1.));
} }
DrumParametersBraid() = delete; //no default initialization allowed DrumParametersBraid() = delete; //no default initialization allowed
DrumParametersBraid(const DrumParametersBraid&) = default; DrumParametersBraid(const DrumParametersBraid&) = default;
DrumParametersBraid& operator=(const DrumParametersBraid&) = default; DrumParametersBraid& operator=(const DrumParametersBraid&) = default;
~DrumParametersBraid() = default; ~DrumParametersBraid() = default;
private: private:
int generate_sublattice(char sublattice, value_t lc) noexcept{ int generate_sublattice(char sublattice, value_t lc) noexcept{
if(sublattice == 'A' or sublattice == 'a'){ if(sublattice == 'A' or sublattice == 'a'){
s0 = Vec2<value_t> (0, -lc); s0 = Vec2<value_t> (0, -lc);
s1 = Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.); s1 = Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.);
s2 = Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.); s2 = Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.);
return 0; return 0;
} }
else if(sublattice == 'B' or sublattice == 'b'){ else if(sublattice == 'B' or sublattice == 'b'){
s0 = -1*(Vec2<value_t> (0, -lc)); s0 = -1*(Vec2<value_t> (0, -lc));
s1 = -1*(Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.)); s1 = -1*(Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.));
s2 = -1*(Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.)); s2 = -1*(Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.));
return 0; return 0;
} }
else else{
return -1; //invalid sublattice identifier return -1; //invalid sublattice identifier
} }
}
public: public:
value_t k0; //onsite prefactor value_t k0; //onsite prefactor
value_t k1; //coupling prefactor 1 value_t k1; //coupling prefactor 1
value_t k2; //coupling prefactor 2 value_t k2; //coupling prefactor 2
value_t c; //damping coefficient value_t c; //damping coefficient
char sublattice; //sublattice the drum is part of (governs electrode geometry) char sublattice; //sublattice the drum is part of (governs electrode geometry)
Vec2<value_t> s0, s1, s2; //vectors connecting to neighbours, ordered according to drum.hpp Vec2<value_t> s0, s1, s2; //vectors connecting to neighbours, ordered according to drum.hpp
Vec2<value_t> position; //position of the drum in real space Vec2<value_t> position; //position of the drum in real space
}; };
template<typename value_t> template<typename value_t>
std::ostream& operator<<(std::ostream& os, const DrumParametersBraid<value_t>& dp) std::ostream& operator<<(std::ostream& os, const DrumParametersBraid<value_t>& dp)
{ {
return os << "Sublattice: " << dp.sublattice << ", k0: " << dp.k0 << ", k1: " << dp.k1 << ", k2: " << dp.k2 << std::endl << "s0: " << dp.s0 << ", s1: " << dp.s1 << ", s2: " << dp.s2; return os << "Sublattice: " << dp.sublattice << ", k0: " << dp.k0 << ", k1: " << dp.k1 << ", k2: " << dp.k2 << std::endl << "s0: " << dp.s0 << ", s1: " << dp.s1 << ", s2: " << dp.s2;
} }
......
projects/braidingTightBinding/include/drum_variables_braid.hpp
View file @ 5f273894
...@@ -8,16 +8,16 @@ ...@@ -8,16 +8,16 @@
//TODO: Generate an interface. It would be trivial, though. //TODO: Generate an interface. It would be trivial, though.
template <typename value_t> template <typename value_t>
class DrumVariablesBraid{ class DrumVariablesBraid{
public: public:
//constructors //constructors
DrumVariablesBraid() = default; DrumVariablesBraid() = default;
DrumVariablesBraid& operator=(const DrumVariablesBraid&) = default; DrumVariablesBraid& operator=(const DrumVariablesBraid&) = default;
~DrumVariablesBraid() = default; ~DrumVariablesBraid() = default;
value_t t0, t1, t2; //t + dt, i.e. baseline plus vortices along bond 0,1,2 value_t t0, t1, t2; //t + dt, i.e. baseline plus vortices along bond 0,1,2
value_t V; //Coupling to central electrode value_t V; //Coupling to central electrode
value_t x, xdot; //Current Elongation (position) and velocity value_t x, xdot; //Current Elongation (position) and velocity
value_t x_temp, xdot_temp; //used as calculation arguments in rk steps value_t x_temp, xdot_temp; //used as calculation arguments in rk steps
}; };
#endif #endif
projects/braidingTightBinding/include/force_braid.hpp
View file @ 5f273894
...@@ -9,35 +9,35 @@ ...@@ -9,35 +9,35 @@
template <typename value_t, typename params_t, typename vars_t, typename buffer_t> template <typename value_t, typename params_t, typename vars_t, typename buffer_t>
class ForceBraid: public Force<value_t, params_t, vars_t, buffer_t>{ class ForceBraid: public Force<value_t, params_t, vars_t, buffer_t>{
public: public:
ForceBraid() = default; ForceBraid() = default;
~ForceBraid() = default; ~ForceBraid() = default;
value_t operator()( value_t operator()(
const Drum<value_t, params_t, vars_t, buffer_t>& drum, //Drum we're calculating the force on const Drum<value_t, params_t, vars_t, buffer_t>& drum, //Drum we're calculating the force on
const Drum<value_t, params_t, vars_t, buffer_t>& neighbour0, //Neighbour 0 const Drum<value_t, params_t, vars_t, buffer_t>& neighbour0, //Neighbour 0
const Drum<value_t, params_t, vars_t, buffer_t>& neighbour1, //Neighbour 1 const Drum<value_t, params_t, vars_t, buffer_t>& neighbour1, //Neighbour 1
const Drum<value_t, params_t, vars_t, buffer_t>& neighbour2, //Neighbour 2 const Drum<value_t, params_t, vars_t, buffer_t>& neighbour2, //Neighbour 2
const value_t time) const noexcept final override const value_t time) const noexcept final override
{ {
//fetch data //fetch data
//note that no copying is done here, these are all references //note that no copying is done here, these are all references
//TODO: implement drive //TODO: implement drive
const params_t& drum_params = drum.get_parameters(); //parameters of drum const params_t& drum_params = drum.get_parameters(); //parameters of drum
const vars_t& drum_vars = drum.get_variables(); //variables of drum const vars_t& drum_vars = drum.get_variables(); //variables of drum
const vars_t& n0_vars = neighbour0.get_variables(); //variables of neighbour 0 const vars_t& n0_vars = neighbour0.get_variables(); //variables of neighbour 0
const vars_t& n1_vars = neighbour1.get_variables(); //variables of neighbour 1 const vars_t& n1_vars = neighbour1.get_variables(); //variables of neighbour 1
const vars_t& n2_vars = neighbour2.get_variables(); //variables of neighbour 2 const vars_t& n2_vars = neighbour2.get_variables(); //variables of neighbour 2
value_t part1 = - drum_params.k0 * drum_vars.x_temp; value_t part1 = - drum_params.k0 * drum_vars.x_temp;
value_t part2 = - drum_params.c * drum_vars.xdot_temp; value_t part2 = - drum_params.c * drum_vars.xdot_temp;
value_t part3 = 0.; value_t part3 = 0.;
part3 -= drum_params.k1 + drum_params.k2 * n0_vars.t0 * n0_vars.x_temp; part3 -= drum_params.k1 + drum_params.k2 * n0_vars.t0 * n0_vars.x_temp;
part3 -= drum_params.k1 + drum_params.k2 * n1_vars.t1 * n1_vars.x_temp; part3 -= drum_params.k1 + drum_params.k2 * n1_vars.t1 * n1_vars.x_temp;
part3 -= drum_params.k1 + drum_params.k2 * n2_vars.t2 * n2_vars.x_temp; part3 -= drum_params.k1 + drum_params.k2 * n2_vars.t2 * n2_vars.x_temp;
return part1 + part2 + part3; return part1 + part2 + part3;
} }
}; };
......
projects/braidingTightBinding/include/grabber.hpp
View file @ 5f273894
...@@ -7,134 +7,137 @@ ...@@ -7,134 +7,137 @@
template<typename value_t> template<typename value_t>
struct DataStorage{ struct DataStorage{
DataStorage(const size_t num_drums) DataStorage(const size_t num_drums)
{ {
x.reserve(num_drums); x.reserve(num_drums);
xdot.reserve(num_drums); xdot.reserve(num_drums);
t0.reserve(num_drums); t0.reserve(num_drums);
t1.reserve(num_drums); t1.reserve(num_drums);
t2.reserve(num_drums); t2.reserve(num_drums);
drive.reserve(num_drums); drive.reserve(num_drums);
} }
~DataStorage() = default; ~DataStorage() = default;
value_t time; value_t time;
std::vector<value_t> x; std::vector<value_t> x;
std::vector<value_t> xdot; std::vector<value_t> xdot;
std::vector<value_t> t0; std::vector<value_t> t0;
std::vector<value_t> t1; std::vector<value_t> t1;
std::vector<value_t> t2; std::vector<value_t> t2;
std::vector<value_t> drive; std::vector<value_t> drive;
}; };
template <typename value_t, typename drum_t> template <typename value_t, typename drum_t>
class Grabber{ class Grabber{
public: public:
Grabber(const size_t grab_every, const std::string params_file, const std::string adjacency_file, const std::string dynamic_file) Grabber(const size_t grab_every, const std::string params_file, const std::string adjacency_file, const std::string dynamic_file)
: grab_every_(grab_every), grab_cnt_(0), par_f_(params_file), adj_f_(adjacency_file), dyn_f_(dynamic_file) {} : grab_every_(grab_every), grab_cnt_(0), par_f_(params_file), adj_f_(adjacency_file), dyn_f_(dynamic_file) {}
Grabber() = delete; //no default initialization Grabber() = delete; //no default initialization
Grabber(const Grabber&) = default; Grabber(const Grabber&) = default;
~Grabber() = default; ~Grabber() = default;
//call in system constructor //call in system constructor
void init(const value_t t_end, const value_t dt, const std::vector<drum_t>& drums, const std::vector<std::vector<int> >& adjacency) noexcept void init(const value_t t_end, const value_t dt, const std::vector<drum_t>& drums, const std::vector<std::vector<int> >& adjacency) noexcept
{ {
drums_ = drums; drums_ = drums;
adjacency_ = adjacency; adjacency_ = adjacency;
t_end_ = t_end; t_end_ = t_end;
dt_ = dt; dt_ = dt;
} }
//call after each simulation step, checks itself if it should save //call after each simulation step, checks itself if it should save
bool grab(const std::vector<drum_t>& drums, const value_t time) bool grab(const std::vector<drum_t>& drums, const value_t time)
{ {
if(grab_cnt_++ % grab_every_ == 0){ if(grab_cnt_++ % grab_every_ == 0){
data_.push_back(DataStorage<value_t>(drums_.size())); data_.push_back(DataStorage<value_t>(drums_.size()));
data_.back().time = time; data_.back().time = time;
for(auto d: drums){ for(auto d: drums){
data_.back().x.push_back(d.get_variables().x); data_.back().x.push_back(d.get_variables().x);
data_.back().xdot.push_back(d.get_variables().xdot); data_.back().xdot.push_back(d.get_variables().xdot);
data_.back().t0.push_back(d.get_variables().t0); data_.back().t0.push_back(d.get_variables().t0);
data_.back().t1.push_back(d.get_variables().t1); data_.back().t1.push_back(d.get_variables().t1);
data_.back().t2.push_back(d.get_variables().t2); data_.back().t2.push_back(d.get_variables().t2);
data_.back().drive.push_back(d.get_variables().V); data_.back().drive.push_back(d.get_variables().V);
} }
return true; return true;
} }
else else{
return false; return false;
} }
}
//prints the data out. call at end of it all. //prints the data out. call at end of it all.
bool save(){ bool save(){
//TODO: Check for file open //TODO: Check for file open
//TODO: Overload operator<< for DrumParameters and DataStorage. //TODO: Overload operator<< for DrumParameters and DataStorage.
//print parameters //print parameters
std::fstream par (par_f_, par.out); std::fstream par (par_f_, par.out);
par << "# The last drum is to be ignored\n"; par << "# The last drum is to be ignored\n";
par << "# k0 \t k1 \t k2 \t c \t pos_x \t pos_y \t sublattice\n"; par << "# k0 \t k1 \t k2 \t c \t pos_x \t pos_y \t sublattice\n";
for(auto d: drums_){ for(auto d: drums_){
par << d.get_parameters().k0 << " \t"; par << d.get_parameters().k0 << " \t";
par << d.get_parameters().k1 << " \t"; par << d.get_parameters().k1 << " \t";
par << d.get_parameters().k2 << " \t"; par << d.get_parameters().k2 << " \t";
par << d.get_parameters().c << " \t"; par << d.get_parameters().c << " \t";
par << d.get_parameters().position.x() << " \t"; par << d.get_parameters().position.x() << " \t";
par << d.get_parameters().position.y() << " \t"; par << d.get_parameters().position.y() << " \t";
par << d.get_parameters().sublattice << "\n"; par << d.get_parameters().sublattice << "\n";
} }
par.close(); par.close();
//print adjacency vectors //print adjacency vectors
std::fstream adj (adj_f_, adj.out); std::fstream adj (adj_f_, adj.out);
adj << "# <-1> means <no neighbour here>\n"; adj << "# <-1> means <no neighbour here>\n";
adj << "# The last drum is to be ignored\n"; adj << "# The last drum is to be ignored\n";
adj << "# n0 \t n1 \t n2\n"; adj << "# n0 \t n1 \t n2\n";
for(auto v: adjacency_){ for(auto v: adjacency_){
for(auto i: v){ for(auto i: v){
if(i == drums_.size()-1) if(i == drums_.size()-1){
adj << -1 << " \t"; //no neighbour present adj << -1 << " \t"; //no neighbour present
else }
adj << i << " \t"; //neighbour else{
} adj << i << " \t"; //neighbour
adj << "\n"; }
} }
adj.close(); adj << "\n";
}
adj.close();
std::fstream dyn (dyn_f_, dyn.out); std::fstream dyn (dyn_f_, dyn.out);
dyn << "# If adjacency of a bond is -1 in " << adj_f_; dyn << "# If adjacency of a bond is -1 in " << adj_f_;
dyn << ", the coupling is to be ignored\n"; dyn << ", the coupling is to be ignored\n";
dyn << "# The last drum is to be ignored\n"; dyn << "# The last drum is to be ignored\n";
dyn << "# Number of drums : " << drums_.size() << "\n"; dyn << "# Number of drums : " << drums_.size() << "\n";
dyn << "# t \t drum_index \t x \t xdot \t t0 \t t1 \t t2 \t drive\n"; dyn << "# t \t drum_index \t x \t xdot \t t0 \t t1 \t t2 \t drive\n";
for(auto entry: data_){ for(auto entry: data_){
for(size_t i = 0; i < entry.x.size(); ++i){ for(size_t i = 0; i < entry.x.size(); ++i){
dyn << entry.time << " \t"; dyn << entry.time << " \t";
dyn << i << " \t"; dyn << i << " \t";
dyn << entry.x[i] << " \t"; dyn << entry.x[i] << " \t";
dyn << entry.xdot[i] << " \t"; dyn << entry.xdot[i] << " \t";
dyn << entry.t0[i] << " \t"; dyn << entry.t0[i] << " \t";
dyn << entry.t1[i] << " \t"; dyn << entry.t1[i] << " \t";
dyn << entry.t2[i] << " \t"; dyn << entry.t2[i] << " \t";
dyn << entry.drive[i] << "\n"; dyn << entry.drive[i] << "\n";
} }
} }
dyn.close(); dyn.close();
return true; return true;
} }
private: private:
size_t grab_every_; //grab data every grab_every_ steps size_t grab_every_; //grab data every grab_every_ steps
size_t grab_cnt_; //grab count size_t grab_cnt_; //grab count
value_t t_end_, dt_; //maximal simulation time and timestep value_t t_end_, dt_; //maximal simulation time and timestep
std::vector<std::vector<int> > adjacency_; //adjacency vector std::vector<std::vector<int> > adjacency_; //adjacency vector
std::vector<drum_t> drums_; //drums in initial configuration, grab params from here std::vector<drum_t> drums_; //drums in initial configuration, grab params from here
std::list<DataStorage<value_t> > data_; //here goes the dynamical data std::list<DataStorage<value_t> > data_; //here goes the dynamical data
//filenames //filenames
std::string par_f_; //parameters file std::string par_f_; //parameters file
std::string adj_f_; //adjacency file std::string adj_f_; //adjacency file
std::string dyn_f_; //dynamics file (this is the gold) std::string dyn_f_; //dynamics file (this is the gold)
}; };
#endif #endif
projects/braidingTightBinding/include/matrix_element_calculator_braid.hpp
View file @ 5f273894
...@@ -4,32 +4,32 @@ ...@@ -4,32 +4,32 @@
template <typename value_t, typename params_t, typename vars_t, typename drum_t> template <typename value_t, typename params_t, typename vars_t, typename drum_t>
class MatrixElementCalculatorBraid: public MatrixElementCalculator<value_t, params_t, vars_t, drum_t>{ class MatrixElementCalculatorBraid: public MatrixElementCalculator<value_t, params_t, vars_t, drum_t>{
public: public:
MatrixElementCalculatorBraid() = default; MatrixElementCalculatorBraid() = default;
~MatrixElementCalculatorBraid() = default; ~MatrixElementCalculatorBraid() = default;
//diagonal element at (index, index) //diagonal element at (index, index)
value_t operator()(const size_t index, const std::vector<drum_t>& drums) const noexcept final override value_t operator()(const size_t index, const std::vector<drum_t>& drums) const noexcept final override
{ {
return drums[index].get_parameters().k0; return drums[index].get_parameters().k0;
} }
//coupling elements at (index1, index2) //coupling elements at (index1, index2)
//for neighbour 0 //for neighbour 0
value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept final override value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept final override
{ {
return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t0; return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t0;
} }
//for neighbour 1 //for neighbour 1
value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int) const noexcept final override value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int) const noexcept final override
{ {
return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t1; return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t1;
} }
//for neighbour 2 //for neighbour 2
value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int, const int) const noexcept final override value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int, const int) const noexcept final override
{ {
return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t2; return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t2;
} }
}; };
#endif #endif
projects/braidingTightBinding/include/vortex.hpp
View file @ 5f273894
...@@ -6,43 +6,43 @@ ...@@ -6,43 +6,43 @@
template <typename value_t> template <typename value_t>
class Vortex{ class Vortex{
public: public:
//constructors //constructors
Vortex(const Vec2<value_t> position, const value_t l0, const value_t alpha, const value_t delta, const value_t a) noexcept: position_(position), l0_(l0), alpha_(alpha), delta_(delta), K_(value_t(4.*3.141592653589793/(3.*std::sqrt(3)*a)), value_t(0.)) {} Vortex(const Vec2<value_t> position, const value_t l0, const value_t alpha, const value_t delta, const value_t a) noexcept: position_(position), l0_(l0), alpha_(alpha), delta_(delta), K_(value_t(4.*3.141592653589793/(3.*std::sqrt(3)*a)), value_t(0.)) {}
Vortex(const Vec2<value_t> position, const Vortex& o) noexcept: position_(position), l0_(o.l0_), alpha_(o.alpha_), delta_(o.delta_), K_(o.K_) {} Vortex(const Vec2<value_t> position, const Vortex& o) noexcept: position_(position), l0_(o.l0_), alpha_(o.alpha_), delta_(o.delta_), K_(o.K_) {}
Vortex() = default; Vortex() = default;
Vortex(const Vortex&) = default; Vortex(const Vortex&) = default;
Vortex& operator=(const Vortex&) = default; Vortex& operator=(const Vortex&) = default;
~Vortex() = default; ~Vortex() = default;
//access //access
value_t l0() const noexcept { return l0_; } value_t l0() const noexcept { return l0_; }
value_t alpha() const noexcept { return alpha_; } value_t alpha() const noexcept { return alpha_; }
value_t delta() const noexcept { return delta_; } value_t delta() const noexcept { return delta_; }
Vec2<value_t>& position() noexcept { return position_; } Vec2<value_t>& position() noexcept { return position_; }
const Vec2<value_t>& position() const noexcept { return position_; } const Vec2<value_t>& position() const noexcept { return position_; }
//modifiers //modifiers
void set_l0(value_t l0) noexcept { l0_ = l0; } void set_l0(value_t l0) noexcept { l0_ = l0; }
void set_alpha(value_t alpha) noexcept { alpha_ = alpha; } void set_alpha(value_t alpha) noexcept { alpha_ = alpha; }
void set_delta(value_t delta) noexcept { delta_ = delta; } void set_delta(value_t delta) noexcept { delta_ = delta; }
void set_position(const Vec2<value_t>& position) noexcept { position_ = position; } void set_position(const Vec2<value_t>& position) noexcept { position_ = position; }
void move_by(const Vec2<value_t>& translation) noexcept { position_ += translation; } void move_by(const Vec2<value_t>& translation) noexcept { position_ += translation; }
void move_to(const Vec2<value_t>& position) noexcept { position_ = position; } void move_to(const Vec2<value_t>& position) noexcept { position_ = position; }
//functional //functional
//there is no divide_by_zero check for efficiency reasons. l0_ can not be zero upon call. //there is no divide_by_zero check for efficiency reasons. l0_ can not be zero upon call.
std::complex<value_t> distortion(const Vec2<value_t>& v) const std::complex<value_t> distortion(const Vec2<value_t>& v) const
{ {
return std::polar<value_t>(delta_*std::tanh(v.r_wrt(position_)/l0_), alpha_ - v.phi_wrt(position_)); return std::polar<value_t>(delta_*std::tanh(v.r_wrt(position_)/l0_), alpha_ - v.phi_wrt(position_));
} }
std::complex<value_t> kekule(const Vec2<value_t>& v, const Vec2<value_t>& s) const{ std::complex<value_t> kekule(const Vec2<value_t>& v, const Vec2<value_t>& s) const{
return std::polar<value_t>(1., K_*(s + 2.*(v - position_))); return std::polar<value_t>(1., K_*(s + 2.*(v - position_)));
} }
private: private:
Vec2<value_t> position_; Vec2<value_t> position_;
const Vec2<value_t> K_; const Vec2<value_t> K_;
value_t l0_, alpha_, delta_; value_t l0_, alpha_, delta_;
}; };
#endif #endif
projects/braidingTightBinding/lib/coupler.hpp
View file @ 5f273894
...@@ -4,14 +4,14 @@ ...@@ -4,14 +4,14 @@
template <typename value_t, typename drum_t> template <typename value_t, typename drum_t>
class Coupler{ class Coupler{
public: public:
//constructors //constructors
Coupler() = default; Coupler() = default;
~Coupler() = default; ~Coupler() = default;
virtual void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept = 0; virtual void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept = 0;
virtual void step(value_t dt) noexcept = 0; virtual void step(value_t dt) noexcept = 0;
virtual value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept = 0; virtual value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept = 0;
}; };
#endif #endif
projects/braidingTightBinding/lib/diagonalizer.hpp
View file @ 5f273894
...@@ -18,74 +18,78 @@ extern "C" void dsyev_( ...@@ -18,74 +18,78 @@ extern "C" void dsyev_(
//TODO: Rule of 7 conformity (nontrivial heap resource) //TODO: Rule of 7 conformity (nontrivial heap resource)
class Diagonalizer{ class Diagonalizer{
public: public:
Diagonalizer() = default; Diagonalizer() = default;
~Diagonalizer() = default; ~Diagonalizer() = default;
std::vector<double> ev(const std::vector<double>& matrix, const size_t N){ std::vector<double> ev(const std::vector<double>& matrix, const size_t N){
//prepare memory //prepare memory
matrix_ = matrix; matrix_ = matrix;
if(eigenvalues_.size() != N){ if(eigenvalues_.size() != N){
eigenvalues_.clear(); eigenvalues_.clear();
eigenvalues_.reserve(N - eigenvalues_.capacity()); eigenvalues_.reserve(N - eigenvalues_.capacity());
} }
for(size_t i = 0; i < N; ++i) for(size_t i = 0; i < N; ++i){
eigenvalues_.push_back(0.); eigenvalues_.push_back(0.);
}
//prepare workspace
N_ = N; //prepare workspace
prepare_workspace('N'); N_ = N;
prepare_workspace('N');
dsyev_('N', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
dsyev_('N', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
if(info_ != 0)
throw("Diagonalization failed!"); if(info_ != 0){
throw("Diagonalization failed!");
//clean up }
delete[] work_;
//clean up
return eigenvalues_; delete[] work_;
}
return eigenvalues_;
std::pair<std::vector<double>, std::vector<double> > evv(const std::vector<double>& matrix, const size_t N){ }
//prepare memory
matrix_ = matrix; std::pair<std::vector<double>, std::vector<double> > evv(const std::vector<double>& matrix, const size_t N){
if(eigenvalues_.size() != N){ //prepare memory
eigenvalues_.clear(); matrix_ = matrix;
eigenvalues_.reserve(N - eigenvalues_.capacity()); if(eigenvalues_.size() != N){
} eigenvalues_.clear();
for(size_t i = 0; i < N; ++i) eigenvalues_.reserve(N - eigenvalues_.capacity());
eigenvalues_.push_back(0.); }
for(size_t i = 0; i < N; ++i){
//prepare workspace eigenvalues_.push_back(0.);
N_ = N; }
prepare_workspace('V');
//prepare workspace
dsyev_('V', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_); N_ = N;
prepare_workspace('V');
if(info_ != 0)
throw("Diagonalization failed!"); dsyev_('V', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
//clean up if(info_ != 0){
delete[] work_; throw("Diagonalization failed!");
}
return std::pair<std::vector<double>, std::vector<double> > (eigenvalues_, matrix_);
} //clean up
delete[] work_;
private:
void prepare_workspace(char type){ return std::pair<std::vector<double>, std::vector<double> > (eigenvalues_, matrix_);
dsyev_(type, 'L', N_, matrix_.data(), N_, eigenvalues_.data(), &dwork_, -1, info_); }
lwork_ = static_cast<int>(dwork_);
work_ = new double[lwork_]; private:
} void prepare_workspace(char type){
dsyev_(type, 'L', N_, matrix_.data(), N_, eigenvalues_.data(), &dwork_, -1, info_);
int info_; lwork_ = static_cast<int>(dwork_);
double dwork_; work_ = new double[lwork_];
int lwork_; }
int N_; //linear matrix dimension
std::vector<double> matrix_; //the matrix int info_;
std::vector<double> eigenvalues_; //eigenvalues double dwork_;
double* work_; int lwork_;
int N_; //linear matrix dimension
std::vector<double> matrix_; //the matrix
std::vector<double> eigenvalues_; //eigenvalues
double* work_;
}; };
#endif #endif