diff --git a/README.md b/README.md index 1b414ee8..afdd7cb5 100644 --- a/README.md +++ b/README.md @@ -1,44 +1,47 @@ # SP4E Homeworks Students: * Lars B. * Bertil T. In this repository, we will keep all homeworks for SP4E ## Homework 1 See folder named *hw1-conjugate-gradient*. #### Requirements * Numpy * Scipy * Matplotlib * Argparse The code as been tested with Python 3. #### Usage Run the *main_minimize.py* with the following arguments: * -A < matrix elements in row major order > * -b < vector elements > * -x0 < initial guess elements > * -method < CG-scipy or CG-ours > * -plot < True or False > #### Example The following command will run the program with the coefficients given in Exercise 1: python3 main_minimize.py -A 8 0 2 6 -b 0 1 -x0 4 4 -method CG-scipy -plot True NB: Note that the quadratic function in Exercise 1 is implemented with a multiplicative factor 1/2 in order to be consistent with Exercise 2. The following command will run the program with a 2x2 s.p.d. matrix A and 2-dim vector b using our self-implemented conjugate gradient method for solving the LSE Ax=b: python3 main_minimize.py -A 3 0 0 4 -b 4 5 -x0 12 12 -plot True -method CG-ours NB: Note that A must be s.p.d. in order for the conjugate gradient method to correctly solve Ax=b. Thus, the matrix A in Exercise 1 should not be used to compare the Scipy version (CG-scipy) againts our version (CG-ours). In the output of Exercise 1, you will see that it does not solve the LSE Ax=b (i.e., the residual outputted is high). NB: Plotting can only be accomplished in 1D (A and b are both scalars) or 2D (A is a 2x2 matrix and b is a 2-dim vector) as it is not possible to plot higher dimensional problems in a way that makes sense. ## Homework 2 -See the homework2/ folder and the instructions within its own README file. +See the **homework2/** folder and the instructions within its own README file. + +## Homework 3 +See the **hw3-heat-fft/** folder and the instructions within its own README file. diff --git a/hw3-heat-fft/CMakeLists.txt b/hw3-heat-fft/CMakeLists.txt index 8d4c0991..ed7ededc 100644 --- a/hw3-heat-fft/CMakeLists.txt +++ b/hw3-heat-fft/CMakeLists.txt @@ -1,81 +1,93 @@ cmake_minimum_required (VERSION 3.1) project (Particles) cmake_policy(VERSION 3.3) set(CMAKE_CXX_STANDARD 14) +option (USE_FFTW "Use FFTW part of the code" ON) + ############################################################### # Include FFTW library and header files ############################################################### -# Not the most elegant way -set(FFTW_LIBRARIES "/usr/lib/x86_64-linux-gnu/libfftw3.so") -set(FFTW_INCLUDES "/usr/include") +if (USE_FFTW) + # Not the most elegant way, but it works! + set(FFTW_LIBRARIES "/usr/lib/x86_64-linux-gnu/libfftw3.so") + set(FFTW_INCLUDES "/usr/include") + + # Alternative way that does not currently work: + #set(FFTW_LIBRARIES CACHE PATH "library where to search libfftw3") + #find_library(FFTW_LIBRARY libfftw3 ${FFTW_LIBRARIES} /usr/lib) + #set(FFTW_INCLUDES CACHE PATH "path where to search fftw include files") + #find_path(FFTW_INCLUDE fftw3.h ${FFTW_INCLUDES} /usr/include/) -# Alternative way that does not currently work: -#set(FFTW_LIBRARIES CACHE PATH "library where to search libfftw3") -#find_library(FFTW_LIBRARY libfftw3 ${FFTW_LIBRARIES} /usr/lib) -#set(FFTW_INCLUDES CACHE PATH "path where to search fftw include files") -#find_path(FFTW_INCLUDE fftw3.h ${FFTW_INCLUDES} /usr/include/) + include_directories("${FFTW_INCLUDES}") -include_directories("${FFTW_INCLUDES}") -message("FFTW_LIBRARIES = ${FFTW_LIBRARIES}") -message("FFTW_INCLUDES = ${FFTW_INCLUDES}") + message("FFTW_LIBRARIES = ${FFTW_LIBRARIES}") + message("FFTW_INCLUDES = ${FFTW_INCLUDES}") +endif(USE_FFTW) ################################################################ # libpart ################################################################ add_library(part compute_boundary.cc compute_verlet_integration.cc particle.cc planet.cc compute_gravity.cc csv_reader.cc particles_factory_interface.cc planets_factory.cc compute_contact.cc compute_kinetic_energy.cc csv_writer.cc system.cc compute_energy.cc compute_potential_energy.cc ping_pong_ball.cc material_point.cc system_evolution.cc ping_pong_balls_factory.cc compute_interaction.cc compute_temperature.cc material_points_factory.cc ) add_executable(particles main.cc) target_link_libraries(particles part) ################################################################ # Google test ################################################################ add_subdirectory(googletest) add_executable(test_kepler test_kepler.cc) -add_executable(test_fft test_fft.cc) + +if (USE_FFTW) + add_executable(test_fft test_fft.cc) +endif(USE_FFTW) + target_link_libraries(test_kepler part gtest_main gtest pthread) -target_link_libraries(test_fft part gtest_main gtest ${FFTW_LIBRARIES} pthread) + +if (USE_FFTW) + target_link_libraries(test_fft part gtest_main gtest ${FFTW_LIBRARIES} pthread) +endif(USE_FFTW) ################################################################ # Doxygen ################################################################ find_package(Doxygen) if (DOXYGEN_FOUND) # to set other options, read: https://cmake.org/cmake/help/v3.9/module/FindDoxygen.html doxygen_add_docs( doxygen ${PROJECT_SOURCE_DIR} COMMENT "Generate html pages" ) add_custom_target(doc DEPENDS doxygen) endif(DOXYGEN_FOUND) diff --git a/hw3-heat-fft/README.md b/hw3-heat-fft/README.md index c028bc33..5c922c90 100644 --- a/hw3-heat-fft/README.md +++ b/hw3-heat-fft/README.md @@ -1 +1,15 @@ # Homework 3 - Heat Equation with FFT + +### Important note: +The option USE_FFTW in CMake must be ON in order to use the FFT part of the code! + +### Comment about the particle code organization: +We now have a different type of "particle" in addition to Planet and PingPongBall. We call this type of particle a **MaterialPoint** (it derives from the Particle.hh class, just like Planet and PingPongBall) + +The MaterialPoint has two properties; **temperature** and **heat rate**, which are both obtainable through their respective get-functions. + +The **MaterialPointsFactory** derives from the general ParticlesFactoryInterface whcih hold the list/vector of particles (in this case MaterialPoints) and the simulation interface. The createSimulation function in MaterialPointsFactory calls **ComputeTemperature** which solves the heat equation on a square (NB!) grid of particles. + +The **Matrix** class (or struct to be precise) is a general C++ class for storing square matrices. + +The **FFT** class (or struct to be precise) contains our implemented wrapping interfaces to Forward DFT, Inverse DFT and for computing FFT frequencies. diff --git a/hw3-heat-fft/fft.hh b/hw3-heat-fft/fft.hh index a1b74179..4bf14268 100644 --- a/hw3-heat-fft/fft.hh +++ b/hw3-heat-fft/fft.hh @@ -1,86 +1,184 @@ #ifndef FFT_HH #define FFT_HH /* ------------------------------------------------------ */ #include "matrix.hh" #include "my_types.hh" #include /* ------------------------------------------------------ */ struct FFT { static Matrix transform(Matrix& m); static Matrix itransform(Matrix& m); static Matrix> computeFrequencies(int size); }; /* ------------------------------------------------------ */ inline Matrix FFT::transform(Matrix& m_in) { int N = m_in.rows(); - - std::cout << "N = " << N << std::endl; + //std::cout << "DEBUG: N = " << N << std::endl; // Declare/initialize in and out for FFTW std::vector in(N*N); std::vector out(N*N); fftw_plan plan = fftw_plan_dft_2d( N, N, reinterpret_cast(&in[0]), reinterpret_cast(&out[0]), FFTW_FORWARD, FFTW_ESTIMATE ); // Fill "in" with input-matrix values for (auto&& entry : index(m_in)) { // Get i and j index of matrix entry int i = std::get<0>(entry); int j = std::get<1>(entry); // Value of matrix entry at (i,j) complex val = std::get<2>(entry); // Create 1D index k from 2D indices (i,j) int k = i + j * N; // Insert value into in in[k] = val; } // Execute FFTW fftw_execute(plan); // Fill output-matrix with "out"-values Matrix m_out(N); for (auto&& entry : index(m_out)) { int i = std::get<0>(entry); int j = std::get<1>(entry); complex& v = std::get<2>(entry); int k = i + j * N; v = out[k]; } // Destroy FFTW fftw_destroy_plan(plan); return m_out; } /* ------------------------------------------------------ */ inline Matrix FFT::itransform(Matrix& m_in) { + int N = m_in.rows(); + + // Declare/initialize in and out for FFTW + std::vector in(N*N); + std::vector out(N*N); + + fftw_plan plan = fftw_plan_dft_2d( N, + N, + reinterpret_cast(&in[0]), + reinterpret_cast(&out[0]), + FFTW_BACKWARD, + FFTW_ESTIMATE ); + + // Fill "in" with input-matrix values + for (auto&& entry : index(m_in)) { + + // Get i and j index of matrix entry + int i = std::get<0>(entry); + int j = std::get<1>(entry); + + // Value of matrix entry at (i,j) + complex val = std::get<2>(entry); + + // Create 1D index k from 2D indices (i,j) + int k = i + j * N; + + // Insert value into in + in[k] = val / (1.0*N*N); + // Note: the divison by N sq follows from the fact that: + // "FFTW computes unnormalized transforms: a transform followed by its + // inverse will result in the original data multiplied by N (or the + // product of the N’s for each dimension, in multi-dimensions)" + // FFTW documentation + } + + // Execute FFTW + fftw_execute(plan); + + // Fill output-matrix with "out"-values + Matrix m_out(N); + for (auto&& entry : index(m_out)) { + int i = std::get<0>(entry); + int j = std::get<1>(entry); + complex& v = std::get<2>(entry); + int k = i + j * N; + v = out[k]; + } + + // Destroy FFTW + fftw_destroy_plan(plan); + + return m_out; + } /* ------------------------------------------------------ */ /* ------------------------------------------------------ */ inline Matrix> FFT::computeFrequencies(int size) { + // This function returns a matrix with entries equal to + // Freq^2 = (Freq_x)^2 + (Freq_y)^2 + // There is no need to take the square root as only the square of k will be used later + + // This vector contains np.fft.fftfreq(size) * size + std::vector Freq_x(size); + + if (size % 2){ // if N is odd + int half_size = (size-1)/2 + 1; + for (int i = 0; i < half_size; i++){ + Freq_x[i] = i; + } + int j = 0; + half_size--; + for (int i = -half_size; i < 0; i++){ + Freq_x[half_size + 1 + j] = i; + j++; + } + + } else{ // N is even + int half_size = size/2; + for (int i = 0; i < half_size; i++){ + Freq_x[i] = i; + } + int j = 0; + for (int i = -half_size; i < 0; i++){ + Freq_x[half_size + j] = i; + j++; + } + } + + // DEBUG: + // for (int i = 0; i < size; i++) + // std::cout << Freq_x[i] << " "; + // std::cout << std::endl; + + // Fill matrix + Matrix> out(size); + for (auto&& entry : index(out)) { + int i = std::get<0>(entry); + int j = std::get<1>(entry); + std::complex& v = std::get<2>(entry); + std::complex Freq_sq = Freq_x[i]*Freq_x[i] + Freq_x[j]*Freq_x[j]; + v = Freq_sq; + } + return out; } #endif // FFT_HH diff --git a/hw3-heat-fft/fft_generate_test_values.py b/hw3-heat-fft/fft_generate_test_values.py new file mode 100644 index 00000000..68c9db99 --- /dev/null +++ b/hw3-heat-fft/fft_generate_test_values.py @@ -0,0 +1,26 @@ +import numpy as np + +N = 20 +vec = np.fft.fftfreq(N) * N +print(vec) +print("\n") + +mat = np.zeros((N,N)) +for i in range(0, N): + for j in range(0, N): + mat[i,j] = vec[i]**2 + vec[j]**2 +print(mat) +print("\n") +np.savetxt("build/fftfreq_test_values_even.txt", mat) + +N = 19 +vec = np.fft.fftfreq(N) * N +print(vec) +print("\n") +mat = np.zeros((N,N)) +for i in range(0, N): + for j in range(0, N): + mat[i,j] = vec[i]**2 + vec[j]**2 +print(mat) +print("\n") +np.savetxt("build/fftfreq_test_values_odd.txt", mat) diff --git a/hw3-heat-fft/test_fft.cc b/hw3-heat-fft/test_fft.cc index 760652d0..cab18dd6 100644 --- a/hw3-heat-fft/test_fft.cc +++ b/hw3-heat-fft/test_fft.cc @@ -1,39 +1,110 @@ #include "my_types.hh" #include "fft.hh" #include +#include /*****************************************************************/ TEST(FFT, transform) { UInt N = 512; Matrix m(N); Real k = 2 * M_PI / N; for (auto&& entry : index(m)) { int i = std::get<0>(entry); int j = std::get<1>(entry); auto& val = std::get<2>(entry); val = cos(k * i); } Matrix res = FFT::transform(m); for (auto&& entry : index(res)) { int i = std::get<0>(entry); int j = std::get<1>(entry); auto& val = std::get<2>(entry); if (std::abs(val) > 1e-10) std::cout << i << "," << j << " = " << val << std::endl; if (i == 1 && j == 0) ASSERT_NEAR(std::abs(val), N * N / 2, 1e-10); else if (i == N - 1 && j == 0) ASSERT_NEAR(std::abs(val), N * N / 2, 1e-10); else ASSERT_NEAR(std::abs(val), 0, 1e-10); } } /*****************************************************************/ TEST(FFT, inverse_transform) { + + UInt N = 512; + Matrix m(N); + Real k = 2 * M_PI / N; + for (auto&& entry : index(m)) { + int i = std::get<0>(entry); + int j = std::get<1>(entry); + auto& val = std::get<2>(entry); + val = cos(k * i); + } + Matrix dft_m = FFT::transform(m); + Matrix invdft_dft_m = FFT::itransform(dft_m); + + for (auto&& entry : index(invdft_dft_m)) { + int i = std::get<0>(entry); + int j = std::get<1>(entry); + auto& val = std::get<2>(entry); + ASSERT_NEAR(val.real(), cos(k*i), 1e-10); + } + } /*****************************************************************/ + +TEST(FFT, compute_frequencies){ + + double tmp; + double val; + UInt N; + + // These tests are based on values produced by numpy.fft.fftfreq + // which are saved in txt files + // See the python script fft_generate_test_values.py + + /////////////// EVEN TEST /////////////////////// + N = 20; + Matrix> even_freqs = FFT::computeFrequencies(N); + std::ifstream in_even("fftfreq_test_values_even.txt"); + if (!in_even) { + std::cout << "Cannot open even file.\n"; + return; + } + for (int j = 0; j < N; j++) { + for (int i = 0; i < N; i++) { + // // DEBUG: + // std::cout << "(i,j) = (" << i << "," << j << ")" << std::endl; + in_even >> tmp; + val = even_freqs(i,j).real(); + ASSERT_NEAR(val, tmp, 1e-10); + } + } + in_even.close(); + + /////////////// ODD TEST /////////////////////// + N = 19; + Matrix> odd_freqs = FFT::computeFrequencies(N); + std::ifstream in_odd("fftfreq_test_values_odd.txt"); + if (!in_odd) { + std::cout << "Cannot open odd file.\n"; + return; + } + for (int j = 0; j < N; j++) { + for (int i = 0; i < N; i++) { + // // DEBUG: + // std::cout << "(i,j) = (" << i << "," << j << ")" << std::endl; + in_odd >> tmp; + val = odd_freqs(i,j).real(); + ASSERT_NEAR(val, tmp, 1e-10); + } + } + in_odd.close(); + +}