From: Alexander Schmidt Date: Sun, 17 Jul 2016 23:46:19 +0000 (+0200) Subject: Basically it's working X-Git-Url: http://git.treefish.org/~alex/seamulator.git/commitdiff_plain/5baf6f709bb7943c95318e709efa907a0f3f1986 Basically it's working --- diff --git a/CMakeLists.txt b/CMakeLists.txt index 63f645f..986edd4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,12 +2,20 @@ cmake_minimum_required(VERSION 2.8) project(seamulator) +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14") + +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}") + find_package(OpenGL REQUIRED) find_package(GLUT REQUIRED) -include_directories(${OPENGL_INCLUDE_DIRS} ${GLUT_INCLUDE_DIRS}) +find_package(fftw3 REQUIRED) + +include_directories(${OPENGL_INCLUDE_DIRS} ${GLUT_INCLUDE_DIRS} + ${FFTW_INCLUDES}) add_executable(seamulator seamulator.cpp sea.cpp watersurface.cpp surfacepoint.cpp seaview.cpp) -target_link_libraries(seamulator ${OPENGL_LIBRARIES} ${GLUT_LIBRARY}) +target_link_libraries(seamulator ${OPENGL_LIBRARIES} ${GLUT_LIBRARY} + ${FFTW_LIBRARIES} ${FFTW_LIBRARIES_THREADS}) diff --git a/Findfftw3.cmake b/Findfftw3.cmake new file mode 100644 index 0000000..00c3401 --- /dev/null +++ b/Findfftw3.cmake @@ -0,0 +1,22 @@ +# - Find FFTW +# Find the native FFTW includes and library +# +# FFTW_INCLUDES - where to find fftw3.h +# FFTW_LIBRARIES - List of libraries when using FFTW. +# FFTW_FOUND - True if FFTW found. + +if (FFTW_INCLUDES) + # Already in cache, be silent + set (FFTW_FIND_QUIETLY TRUE) +endif (FFTW_INCLUDES) + +find_path (FFTW_INCLUDES fftw3.h) + +find_library (FFTW_LIBRARIES NAMES fftw3) + +# handle the QUIETLY and REQUIRED arguments and set FFTW_FOUND to TRUE if +# all listed variables are TRUE +include (FindPackageHandleStandardArgs) +find_package_handle_standard_args (FFTW DEFAULT_MSG FFTW_LIBRARIES FFTW_INCLUDES) + +mark_as_advanced (FFTW_LIBRARIES FFTW_INCLUDES) diff --git a/sea.cpp b/sea.cpp index 30593f2..b76a9dd 100644 --- a/sea.cpp +++ b/sea.cpp @@ -6,8 +6,8 @@ #include "watersurface.h" -const double Sea::PHILLIPS_CONSTANT{1}; -const double Sea::GRAVITATIONAL_CONSTANT{6.67408e-11}; +const double Sea::PHILLIPS_CONSTANT{0.0000003}; +const double Sea::GRAVITATIONAL_CONSTANT{9.8}; Sea::Sea(WaterSurfacePtr surface) : m_surface{surface}, @@ -18,14 +18,68 @@ Sea::Sea(WaterSurfacePtr surface) : { m_fourierAmplitudes.resize(pow(m_surface->size(), 2)); generateFourierAmplitudes(); + + fftw_init_threads(); + fftw_plan_with_nthreads(4); + + m_fftwIn = (fftw_complex*) + fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2)); + m_fftwOut = (fftw_complex*) + fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2)); + m_fftwPlan = fftw_plan_dft_2d + (m_surface->size(), m_surface->size(), m_fftwIn, m_fftwOut, + FFTW_BACKWARD, FFTW_MEASURE); + + m_startTime = std::chrono::system_clock::now(); +} + +Sea::~Sea() +{ + fftw_destroy_plan(m_fftwPlan); + fftw_free(m_fftwIn); fftw_free(m_fftwOut); +} + +double Sea::getRuntime() const +{ + auto timeNow = std::chrono::system_clock::now(); + auto durationMs = + std::chrono::duration_cast(timeNow - m_startTime); + + return durationMs.count() / 1000.0; } void Sea::update() { + using namespace std::complex_literals; + + const double runtime = getRuntime(); + + for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) { + const int positiveM = (m + m_surface->size()) % m_surface->size(); + + for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) { + const double k = sqrt(pow(spatialFrequencyForIndex(n), 2) + + pow(spatialFrequencyForIndex(m), 2)); + const double omega = sqrt(GRAVITATIONAL_CONSTANT * k); + + std::complex amplitude = + fourierAmplitudeAt(n, m).first * exp(1i * omega * runtime) + + fourierAmplitudeAt(n, m).second * exp(-1i * omega * runtime); + + const int positiveN = (n + m_surface->size()) % m_surface->size(); + int fftwIndex = positiveM + positiveN * m_surface->size(); + + m_fftwIn[fftwIndex][0] = std::real(amplitude); + m_fftwIn[fftwIndex][1] = std::imag(amplitude); + } + } + + fftw_execute(m_fftwPlan); + for (int y = 0; y < m_surface->size(); ++y) { for (int x = 0; x < m_surface->size(); ++x) { m_surface->at(x, y) - .setHeight(((double)std::rand()/(double)RAND_MAX)); + .setHeight(m_fftwOut[y + x * m_surface->size()][0]); } } } @@ -72,4 +126,6 @@ void Sea::generateFourierAmplitudes() spatialFrequencyForIndex(m)); } } + + fourierAmplitudeAt(0, 0) = {0, 0}; } diff --git a/sea.h b/sea.h index 6e24ac3..20f5d63 100644 --- a/sea.h +++ b/sea.h @@ -1,8 +1,11 @@ #pragma once +#include #include #include +#include + #include "seafwd.h" #include "complexpair.h" @@ -12,6 +15,9 @@ class Sea { public: Sea(WaterSurfacePtr surface); + ~Sea(); + Sea(const Sea&) = delete; + Sea& operator=(const Sea&) = delete; void update(); private: @@ -25,10 +31,14 @@ class Sea std::mt19937 m_randomGenerator; std::normal_distribution<> m_normalDistribution; std::vector m_fourierAmplitudes; + fftw_complex *m_fftwIn, *m_fftwOut; + fftw_plan m_fftwPlan; + std::chrono::time_point m_startTime; double phillipsSpectrum(double k_x, double k_y) const; ComplexPair generateFourierAmplitude(double k_x, double k_y); ComplexPair& fourierAmplitudeAt(int n, int m); void generateFourierAmplitudes(); double spatialFrequencyForIndex(int n) const; + double getRuntime() const; }; diff --git a/seamulator.cpp b/seamulator.cpp index 1c1ff67..c5f7936 100644 --- a/seamulator.cpp +++ b/seamulator.cpp @@ -8,12 +8,12 @@ #include "seaview.h" #include "watersurface.h" -const int LATTICE_SIZE{10}; +const int LATTICE_SIZE{128}; const double LATTICE_EXTEND{10}; -const int INIT_WINDOW_POS_X{100}; -const int INIT_WINDOW_POS_Y{100}; -const int INIT_WINDOW_WIDTH{300}; -const int INIT_WINDOW_HEIGHT{300}; +const int INIT_WINDOW_POS_X{50}; +const int INIT_WINDOW_POS_Y{50}; +const int INIT_WINDOW_WIDTH{800}; +const int INIT_WINDOW_HEIGHT{600}; const double INIT_VIEW_DISTANCE{LATTICE_EXTEND * 1.5}; const double INIT_VIEW_AZIMUTH{0}; const double INIT_VIEW_ALTITUDE{M_PI / 4};