diff --git a/PlotScript/Reserve.txt b/PlotScript/Reserve.txt deleted file mode 100644 index bd36799..0000000 --- a/PlotScript/Reserve.txt +++ /dev/null @@ -1,11 +0,0 @@ -4.0 -4.0 -4.0 -0.0 -1.0 -0.0 -1.0 -0.0 -1.0 -4.0 -0.005 diff --git a/PlotScript/Source.txt b/PlotScript/Source.txt deleted file mode 100644 index dbb3d9d..0000000 --- a/PlotScript/Source.txt +++ /dev/null @@ -1,11 +0,0 @@ -32.0 -32.0 -32.0 -0.0 -1.0 -0.0 -1.0 -0.0 -1.0 -32.0 -0.005 diff --git a/PlotScript/convergence.py b/PlotScript/convergence.py deleted file mode 100644 index 6e44551..0000000 --- a/PlotScript/convergence.py +++ /dev/null @@ -1,98 +0,0 @@ -import matplotlib.pyplot as plt -import subprocess - -components = {1: "Ex", 2: "Ey", 3: "Ez", 4: "Bx", 5: "By", 6: "Bz"} -nums_com = {"Ex": 0, "Ey": 1, "Ez": 2, "Bx": 3, "By": 4, "Bz": 5} -input_list = ["Ni", "Nj", "Nk" "ax", "bx", "ay", "by", "az", "bz", "dt", "t"] - - -def execute_cpp(field_1, field_2, field_to_plot, shift_flag): - num_field_1 = nums_com[field_1] - num_field_2 = nums_com[field_2] - num_field_to_plot = nums_com[field_to_plot] - - print("\n" + field_to_plot + ":\n") - - cpp_executable = "src/Release/sample.exe" - args = [cpp_executable, str(num_field_1), str(num_field_2), str(num_field_to_plot), shift_flag] - try: - completed_process = subprocess.run(args, check=True, stdout=subprocess.PIPE) - output = completed_process.stdout.decode("utf-8").strip() - print(output) - return output - except subprocess.CalledProcessError: - print("Error when starting a C++ project") - except FileNotFoundError: - print("sample.exe not found") - - -def update_sources(): - print("Update parameters? \n \ - 1 - Yes \n \ - 2 - No") - select_update = int(input("Number: ")) * (-1) + 2 - if not (select_update == 0 or select_update == 1): - print("Invalid input") - exit(1) - elif select_update == 1: - with open("Source.txt", "w") as file: - for component in input_list: - file.write(input(component + ": ") + "\n") - - -def save_source_into_reserve(): - with open('Source.txt', 'r') as file1, open('Reserve.txt', 'w') as file2: - numbers = [float(line.strip()) for line in file1] - for num in numbers: - file2.write(str(num) + "\n") - - return numbers - - -def reload_source(): - with open('Source.txt', 'w') as file2, open('Reserve.txt', 'r') as file1: - numbers = [float(line.strip()) for line in file1] - for num in numbers: - file2.write(str(num) + "\n") - - -def analyze_convergence(numbers, mult_2, iterations, shifts): - convergences = [] - nums = [] - flag = str(int(shifts)) - for n in range(0, iterations): - with open("Source.txt", "w") as file: - mult_1 = 2 - for num in range(0, 3): - tmp = numbers[num] * (mult_1 ** n) - file.write(str(tmp) + "\n") - - for num in range(3, 9): - file.write(str(numbers[num]) + "\n") - - tmp_n_6 = numbers[9] * (mult_2 ** n) - file.write(str(tmp_n_6) + "\n") - file.write(str(numbers[10]) + "\n") - - convergences.append(float(execute_cpp("Ex", "By", "Ex", flag))) - - convers = [] - for n in range(0, iterations - 1): - convers.append(convergences[n] / convergences[n + 1]) - nums.append(n) - return nums, convers - - -if __name__ == '__main__': - update_sources() - loaded_numbers = save_source_into_reserve() - nums, convergence = analyze_convergence(loaded_numbers, 2, 4, False) - reload_source() - - print(convergence) - plt.plot(nums, convergence) - plt.xlabel('n') - plt.ylabel('error reduction factor') - plt.title("Plot") - plt.grid(True) - plt.show() diff --git a/PlotScript/main.py b/PlotScript/main.py deleted file mode 100644 index bfb536e..0000000 --- a/PlotScript/main.py +++ /dev/null @@ -1,100 +0,0 @@ -import matplotlib.pyplot as plt -import csv -import subprocess - - -components = {1: 'Ex', 2: 'Ey', 3: 'Ez', 4: 'Bx', 5: 'By', 6: 'Bz'} -nums_com = {'Ex': 0, 'Ey': 1, 'Ez': 2, 'Bx': 3, 'By': 4, 'Bz': 5} -shift_flag = "1" - - -def get_plot(field, data_path, size_axis, data_block_height, axis_boundaries, time): - field_num = nums_com[field] + 1 - a = 0. - A = [] - V = [] - dA = (float(axis_boundaries[1]) - float(axis_boundaries[0])) / float(size_axis) - - for n in range(size_axis): - A.append(a) - a += dA - with open(data_path, 'r') as datafile: - plotting = list(csv.reader(datafile, delimiter=';')) - V.extend([float(value) for value in plotting[(field_num - 1) * (data_block_height + 2)]]) - - plt.plot(A, V) - plt.xlabel('Coordinate') - plt.ylabel(components[field_num]) - plt.title(f'Plot {components[field_num]} (Time: {str(time)})') - plt.grid(True) - plt.savefig(f'Plots/plt_{components[field_num]}') - plt.show() - - -def select_parameters(field_E, field_B, source_nums): - num_field_E = nums_com[field_E] - num_field_B = nums_com[field_B] - - if num_field_E == 0 and num_field_B == 5 or num_field_E == 2 and num_field_B == 3: - boundaries = (float(source_nums[3]), float(source_nums[4])) - data_size = source_nums[1] * source_nums[2] - axis_size = source_nums[0] - elif num_field_E == 1 and num_field_B == 5 or num_field_E == 2 and num_field_B == 4: - boundaries = (float(source_nums[5]), float(source_nums[6])) - data_size = source_nums[0] * source_nums[2] - axis_size = source_nums[1] - elif num_field_E == 0 and num_field_B == 4 or num_field_E == 1 and num_field_B == 3: - boundaries = (float(source_nums[7]), float(source_nums[8])) - data_size = source_nums[0] * source_nums[1] - axis_size = source_nums[2] - else: - print('Invalid selected fields') - exit(1) - - return int(axis_size), boundaries, int(data_size) - - -def execute_cpp(field_E, field_B, field_to_plot): - num_field_E = nums_com[field_E] - num_field_B = nums_com[field_B] - num_field_to_plot = nums_com[field_to_plot] - - print("\n" + field_to_plot + ":\n") - - cpp_executable = 'src/Release/sample.exe' - args = [cpp_executable, str(num_field_E), str(num_field_B), str(num_field_to_plot), shift_flag] - try: - subprocess.run(args, check=True) - except subprocess.CalledProcessError: - print('Error when starting a C++ project') - except FileNotFoundError: - print('sample.exe not found') - - -if __name__ == '__main__': - input_list = ['Ni', 'Nj', 'Nk' 'ax', 'bx', 'ay', 'by', 'az', 'bz', 'Nt', 't'] - - print('Update parameters? \n \ - 1 - Yes \n \ - 2 - No') - select_update = int(input('Number: ')) * (-1) + 2 - if not (select_update == 0 or select_update == 1): - print('Invalid input') - exit(1) - - if (select_update): - with open('Source.txt', 'w') as file: - for component in input_list: - file.write(input(component + ': ') + '\n') - - with open('Source.txt', 'r') as file: - numbers = [float(line.strip()) for line in file] - - field_E = 'Ez' - field_B = 'Bx' - field_to_plot = 'Ez' - - execute_cpp(field_E, field_B, field_to_plot) - size_axis, axis_boundaries, block_size = select_parameters(field_E, field_B, numbers) - file = 'OutFile.csv' - get_plot(field_to_plot, file, size_axis, block_size, axis_boundaries, numbers[10]) diff --git a/PlotScript/src/Release/sample.exe b/PlotScript/src/Release/sample.exe deleted file mode 100644 index 2e15bb8..0000000 Binary files a/PlotScript/src/Release/sample.exe and /dev/null differ diff --git a/PlotScript/visualization.py b/PlotScript/visualization.py index f245f16..487c039 100644 --- a/PlotScript/visualization.py +++ b/PlotScript/visualization.py @@ -1,5 +1,8 @@ +import os import glob import subprocess +import platform +import argparse import pandas as pd import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation @@ -9,21 +12,21 @@ def get_animation(component): folder_path = f'OutFiles_{components_num[component]}' file_names_base = sorted(glob.glob(folder_path + '/*.csv')) - file_names = sorted(file_names_base, key=lambda x: int(x.split('\\')[-1].split('.')[0])) + file_names = sorted(file_names_base, key=lambda x: int(x.split(os.path.sep)[-1].split('.')[0])) fig, ax = plt.subplots() def update(frame): - mut = 0 ax.clear() data = pd.read_csv(file_names[frame], sep=';') data = data.apply(lambda x: x.str.replace(',', '.').astype(float) if x.dtype == 'object' else x) ax.imshow(data, cmap='viridis', vmin=-0.07, vmax=0.07, interpolation='none') # RdBu ax.set_title('Frame {}'.format(frame + 1)) - + + if not os.path.exists('animations'): + os.makedirs('animations') + ani = FuncAnimation(fig, update, frames=len(file_names), interval=50) - - plt.show() ani.save(f'animations/animation_{component}.gif', writer='imagemagick') @@ -32,22 +35,54 @@ def get_heatmap(component, iteration): data = data.applymap(lambda x: float(x.replace(',', '.')) if isinstance(x, str) else x) - plt.imshow(data, cmap='viridis', interpolation='nearest') + plt.imshow(data, cmap='viridis', vmin=-0.07, vmax=0.07, interpolation='none') plt.colorbar() - plt.show() + if not os.path.exists('heatmap'): + os.makedirs('heatmap') + + heatmap_path = os.path.join('heatmap', f'heatmap_iteration_{iteration}_component_{component}.png') + plt.savefig(heatmap_path) + + +def execute_cpp(grid_size, iters_num): + + system = platform.system() -def execute_cpp(grid_size, iters_num, single_iteration_flag=True): - cpp_executable = 'src/Release/sample.exe' - args = [cpp_executable, str(grid_size), str(iters_num), str(int(single_iteration_flag))] + if system == 'Windows': + cpp_executable = 'src/Release/sample.exe' + else: + cpp_executable = 'src/Release/sample' + + args = [cpp_executable, str(grid_size), str(iters_num), str(int(True))] try: subprocess.run(args, check=True) except subprocess.CalledProcessError: print('Error when starting a C++ project') except FileNotFoundError: - print('sample.exe not found') + print(f'{cpp_executable} not found') if __name__ == '__main__': - execute_cpp(75, 410, False) - get_animation("Ex") + parser = argparse.ArgumentParser() + + parser.add_argument('--run_cpp', action='store_true', help="Run C++ application before generating data") + parser.add_argument('--function', choices=['heatmap', 'animation'], help="Select function to execute") + parser.add_argument('--iteration', type=int, default=0, help="Iteration number for heatmap") + parser.add_argument('component', type=str, help="Component for analysis") + + group = parser.add_argument_group('conditional arguments') + group.add_argument('--grid_size', type=int, help="Grid size") + group.add_argument('--iters_num', type=int, help="Number of iterations") + + args = parser.parse_args() + + if args.run_cpp: + if not (args.grid_size and args.iters_num): + parser.error("The --run_cpp flag requires --grid_size and --iters_num arguments") + execute_cpp(args.grid_size, args.iters_num) + + if args.function == 'heatmap': + get_heatmap(args.component, args.iteration) + elif args.function == 'animation': + get_animation(args.component) diff --git a/include/Structures.h b/include/Structures.h index 28f21d4..5ac42bc 100644 --- a/include/Structures.h +++ b/include/Structures.h @@ -9,6 +9,7 @@ namespace FDTDconst const double EPS0 = 1.0; const double MU0 = EPS0; const double N = 4.0; + const double PI = 3.14159265358; } namespace FDTDstruct diff --git a/sln/CMakeLists.txt b/sln/CMakeLists.txt index 476847a..291e5c1 100644 --- a/sln/CMakeLists.txt +++ b/sln/CMakeLists.txt @@ -4,6 +4,13 @@ project(FDTD) include_directories(../include ../gtest ../test) +find_package(OpenMP REQUIRED) + +if(OPENMP_FOUND) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}") +endif() # BUILD add_subdirectory(FDTD) diff --git a/sln/sample/CMakeLists.txt b/sln/sample/CMakeLists.txt index 544d079..1dd577d 100644 --- a/sln/sample/CMakeLists.txt +++ b/sln/sample/CMakeLists.txt @@ -2,5 +2,15 @@ file(GLOB hdrs "*.h*" "../../include/*.h") file(GLOB srcs "*.cpp" "../../src/*.cpp" "../../samples/*.cpp") add_executable(sample ${srcs} ${hdrs}) -set_target_properties(sample PROPERTIES RUNTIME_OUTPUT_DIRECTORY "../../PlotScript/src") +if(WIN32) + set(output_dir "../../PlotScript/src") +else() + set(output_dir "../../PlotScript/src/Release") +endif() + +set_target_properties(sample PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${output_dir}) + +if(NOT EXISTS ${output_dir}) + file(MAKE_DIRECTORY ${output_dir}) +endif() diff --git a/src/FDTD.cpp b/src/FDTD.cpp index f763286..1e5125f 100644 --- a/src/FDTD.cpp +++ b/src/FDTD.cpp @@ -1,6 +1,8 @@ #include "FDTD.h" #include +#include +#include //#define __DEBUG_PML__ //#define __OUTPUT_SIGMA__ @@ -55,10 +57,10 @@ void FDTD::set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2], for (int k = bounds_k[0]; k < bounds_k[1]; k++) { EsigmaX(i, j, k) = SGm * - pow((static_cast(dist(i, j, k))) / + std::pow((static_cast(dist(i, j, k))) / static_cast(pml_size_i), FDTDconst::N); BsigmaX(i, j, k) = SGm * - pow((static_cast(dist(i, j, k))) / + std::pow((static_cast(dist(i, j, k))) / static_cast(pml_size_i), FDTDconst::N); } } @@ -75,10 +77,10 @@ void FDTD::set_sigma_y(int bounds_i[2], int bounds_j[2], int bounds_k[2], for (int k = bounds_k[0]; k < bounds_k[1]; k++) { EsigmaY(i, j, k) = SGm * - pow((static_cast(dist(i, j, k))) / + std::pow((static_cast(dist(i, j, k))) / static_cast(pml_size_j), FDTDconst::N); BsigmaY(i, j, k) = SGm * - pow((static_cast(dist(i, j, k))) / + std::pow((static_cast(dist(i, j, k))) / static_cast(pml_size_j), FDTDconst::N); } } @@ -95,10 +97,10 @@ void FDTD::set_sigma_z(int bounds_i[2], int bounds_j[2], int bounds_k[2], for (int k = bounds_k[0]; k < bounds_k[1]; k++) { EsigmaZ(i, j, k) = SGm * - pow((static_cast(dist(i, j, k))) / + std::pow((static_cast(dist(i, j, k))) / static_cast(pml_size_k), FDTDconst::N); BsigmaZ(i, j, k) = SGm * - pow((static_cast(dist(i, j, k))) / + std::pow((static_cast(dist(i, j, k))) / static_cast(pml_size_k), FDTDconst::N); } } @@ -113,7 +115,7 @@ FDTD::FDTD(Parameters _parameters, double _dt, double _pml_percent) : parameters.Nk <= 0 || dt <= 0) { - throw std::exception("ERROR: invalid parameters"); + throw std::invalid_argument("ERROR: invalid parameters"); } Ex = Ey = Ez = Bx = By = Bz = Field(parameters.Ni, parameters.Nj, parameters.Nk); @@ -145,7 +147,7 @@ Field& FDTD::get_field(Component this_field) case Component::BZ: return Bz; - default: throw std::exception("ERROR: Invalid field component"); + default: throw std::logic_error("ERROR: Invalid field component"); } } @@ -159,7 +161,7 @@ std::vector& FDTD::get_current(Component this_current) case Component::JZ: return Jz; - default: throw std::exception("ERROR: Invalid current component"); + default: throw std::logic_error("ERROR: Invalid current component"); } } @@ -176,13 +178,13 @@ void FDTD::update_E(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t) { for (int k = bounds_k[0]; k < bounds_k[1]; k++) { - Ex(i, j, k) = Ex(i, j, k) - 4.0 * M_PI * dt * Jx[t](i, j, k) + + Ex(i, j, k) = Ex(i, j, k) - 4.0 * FDTDconst::PI * dt * Jx[t](i, j, k) + FDTDconst::C * dt * ((Bz(i, j, k) - Bz(i, j - 1, k)) / dy - (By(i, j, k) - By(i, j, k - 1)) / dz); - Ey(i, j, k) = Ey(i, j, k) - 4.0 * M_PI * dt * Jy[t](i, j, k) + + Ey(i, j, k) = Ey(i, j, k) - 4.0 * FDTDconst::PI * dt * Jy[t](i, j, k) + FDTDconst::C * dt * ((Bx(i, j, k) - Bx(i, j, k - 1)) / dz - (Bz(i, j, k) - Bz(i - 1, j, k)) / dx); - Ez(i, j, k) = Ez(i, j, k) - 4.0 * M_PI * dt * Jz[t](i, j, k) + + Ez(i, j, k) = Ez(i, j, k) - 4.0 * FDTDconst::PI * dt * Jz[t](i, j, k) + FDTDconst::C * dt * ((By(i, j, k) - By(i - 1, j, k)) / dx - (Bx(i, j, k) - Bx(i, j - 1, k)) / dy); } @@ -219,7 +221,7 @@ void FDTD::update_B(int bounds_i[2], int bounds_j[2], int bounds_k[2]) double FDTD::PMLcoef(double sigma) { - return exp(-sigma * dt * FDTDconst::C); + return std::exp(-sigma * dt * FDTDconst::C); } void FDTD::update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) @@ -341,24 +343,24 @@ double FDTD::calc_reflection(Field E_start[3], Field B_start[3], { for (int k = pml_size_k; k < parameters.Nk - pml_size_k; k++) { - energy_start += pow(E_start[0](i, j, k), 2.0) + - pow(E_start[1](i, j, k), 2.0) + - pow(E_start[2](i, j, k), 2.0) + - pow(B_start[0](i, j, k), 2.0) + - pow(B_start[1](i, j, k), 2.0) + - pow(B_start[2](i, j, k), 2.0); - - energy_final += pow(E_final[0](i, j, k), 2.0) + - pow(E_final[1](i, j, k), 2.0) + - pow(E_final[2](i, j, k), 2.0) + - pow(B_final[0](i, j, k), 2.0) + - pow(B_final[1](i, j, k), 2.0) + - pow(B_final[2](i, j, k), 2.0); + energy_start += std::pow(E_start[0](i, j, k), 2.0) + + std::pow(E_start[1](i, j, k), 2.0) + + std::pow(E_start[2](i, j, k), 2.0) + + std::pow(B_start[0](i, j, k), 2.0) + + std::pow(B_start[1](i, j, k), 2.0) + + std::pow(B_start[2](i, j, k), 2.0); + + energy_final += std::pow(E_final[0](i, j, k), 2.0) + + std::pow(E_final[1](i, j, k), 2.0) + + std::pow(E_final[2](i, j, k), 2.0) + + std::pow(B_final[0](i, j, k), 2.0) + + std::pow(B_final[1](i, j, k), 2.0) + + std::pow(B_final[2](i, j, k), 2.0); } } } - return sqrt(energy_final / energy_start); + return std::sqrt(energy_final / energy_start); } double FDTD::calc_max_value(Field& field) @@ -372,7 +374,7 @@ double FDTD::calc_max_value(Field& field) { for (int k = 0; k < parameters.Nk - pml_size_k; k++) { - if (fabs(field(i, j, k)) > fabs(max_val)) + if (std::abs(field(i, j, k)) > std::abs(max_val)) { max_val = field(i, j, k); } @@ -386,7 +388,7 @@ std::vector> FDTD::update_fields(const int time) { if (time < 0) { - throw std::exception("ERROR: Invalid update field argument"); + throw std::invalid_argument("ERROR: Invalid update field argument"); } std::vector> return_data; @@ -459,11 +461,11 @@ std::vector> FDTD::update_fields(const int time) }; // Calculation of maximum permittivity - double SGm_x = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) + double SGm_x = -(FDTDconst::N + 1.0) / 2.0 * std::log(FDTDconst::R) / (static_cast(pml_size_i) * parameters.dx); - double SGm_y = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) + double SGm_y = -(FDTDconst::N + 1.0) / 2.0 * std::log(FDTDconst::R) / (static_cast(pml_size_j) * parameters.dy); - double SGm_z = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) + double SGm_z = -(FDTDconst::N + 1.0) / 2.0 * std::log(FDTDconst::R) / (static_cast(pml_size_k) * parameters.dz); // Calculation of permittivity in the cells diff --git a/src/test_FDTD.cpp b/src/test_FDTD.cpp index ee10366..fe0cfea 100644 --- a/src/test_FDTD.cpp +++ b/src/test_FDTD.cpp @@ -126,7 +126,7 @@ void Test_FDTD::set_sign(Component field_E, Component field_B) { sign = 1.0; } - else throw std::exception("ERROR: invalid selected fields"); + else throw std::logic_error("ERROR: invalid selected fields"); } void Test_FDTD::set_axis(Component field_E, Component field_B) { @@ -145,7 +145,7 @@ void Test_FDTD::set_axis(Component field_E, Component field_B) { axis = Axis::Z; } - else throw std::exception("ERROR: invalid selected fields"); + else throw std::logic_error("ERROR: invalid selected fields"); } double Test_FDTD::get_shift(Component _field, double step) {