diff --git a/src/simulation/Simulation.cpp b/src/simulation/Simulation.cpp index 38efca797..13c3e8c82 100644 --- a/src/simulation/Simulation.cpp +++ b/src/simulation/Simulation.cpp @@ -3806,10 +3806,10 @@ void Simulation::BeforeSim() grav->gravity_update_async(); //Get updated buffer pointers for gravity - gravx = grav->gravx; - gravy = grav->gravy; - gravp = grav->gravp; - gravmap = grav->gravmap; + gravx = &grav->gravx[0]; + gravy = &grav->gravy[0]; + gravp = &grav->gravp[0]; + gravmap = &grav->gravmap[0]; } if(emp_decor>0) emp_decor -= emp_decor/25+2; @@ -4045,10 +4045,10 @@ Simulation::Simulation(): //Give air sim references to our data grav->bmap = bmap; //Gravity sim gives us maps to use - gravx = grav->gravx; - gravy = grav->gravy; - gravp = grav->gravp; - gravmap = grav->gravmap; + gravx = &grav->gravx[0]; + gravy = &grav->gravy[0]; + gravp = &grav->gravp[0]; + gravmap = &grav->gravmap[0]; //Create and attach air simulation air = new Air(*this); diff --git a/src/simulation/gravity/Common.cpp b/src/simulation/gravity/Common.cpp index 468a23639..d73f12783 100644 --- a/src/simulation/gravity/Common.cpp +++ b/src/simulation/gravity/Common.cpp @@ -9,44 +9,30 @@ Gravity::Gravity(CtorTag) { - // Allocate full size Gravmaps - unsigned int size = NCELL; - th_ogravmap = new float[size]; - th_gravmap = new float[size]; - th_gravy = new float[size]; - th_gravx = new float[size]; - th_gravp = new float[size]; - gravmap = new float[size]; - gravy = new float[size]; - gravx = new float[size]; - gravp = new float[size]; - gravmask = new unsigned[size]; + th_ogravmap.resize(NCELL); + th_gravmap.resize(NCELL); + th_gravy.resize(NCELL); + th_gravx.resize(NCELL); + th_gravp.resize(NCELL); + gravmap.resize(NCELL); + gravy.resize(NCELL); + gravx.resize(NCELL); + gravp.resize(NCELL); + gravmask.resize(NCELL); } Gravity::~Gravity() { stop_grav_async(); - - delete[] th_ogravmap; - delete[] th_gravmap; - delete[] th_gravy; - delete[] th_gravx; - delete[] th_gravp; - delete[] gravmap; - delete[] gravy; - delete[] gravx; - delete[] gravp; - delete[] gravmask; } void Gravity::Clear() { - int size = NCELL; - std::fill(gravy, gravy + size, 0.0f); - std::fill(gravx, gravx + size, 0.0f); - std::fill(gravp, gravp + size, 0.0f); - std::fill(gravmap, gravmap + size, 0.0f); - std::fill(gravmask, gravmask + size, 0xFFFFFFFF); + std::fill(&gravy[0], &gravy[0] + NCELL, 0.0f); + std::fill(&gravx[0], &gravx[0] + NCELL, 0.0f); + std::fill(&gravp[0], &gravp[0] + NCELL, 0.0f); + std::fill(&gravmap[0], &gravmap[0] + NCELL, 0.0f); + std::fill(&gravmask[0], &gravmask[0] + NCELL, UINT32_C(0xFFFFFFFF)); ignoreNextResult = true; } @@ -86,8 +72,8 @@ void Gravity::gravity_update_async() gravcv.notify_one(); } unsigned int size = NCELL; - membwand(gravy, gravmask, size * sizeof(float), size * sizeof(unsigned)); - membwand(gravx, gravmask, size * sizeof(float), size * sizeof(unsigned)); + membwand(&gravy[0], &gravmask[0], size * sizeof(float), size * sizeof(uint32_t)); + membwand(&gravx[0], &gravmask[0], size * sizeof(float), size * sizeof(uint32_t)); std::fill(&gravmap[0], &gravmap[0] + size, 0.0f); } @@ -244,8 +230,6 @@ void Gravity::gravity_mask() unsigned maskvalue; mask_el *t_mask_el = nullptr; mask_el *c_mask_el = nullptr; - if (!gravmask) - return; memset(checkmap, 0, sizeof(checkmap)); for (int x = 0; x < XCELLS; x++) { diff --git a/src/simulation/gravity/Fft.cpp b/src/simulation/gravity/Fft.cpp index e87ad623c..4d8cac5f8 100644 --- a/src/simulation/gravity/Fft.cpp +++ b/src/simulation/gravity/Fft.cpp @@ -2,19 +2,37 @@ #include "Misc.h" #include #include +#include #include +constexpr auto xblock2 = XCELLS * 2; +constexpr auto yblock2 = YCELLS * 2; +constexpr auto fft_tsize = (xblock2 / 2 + 1) * yblock2; +//NCELL*4 is size of data array, scaling needed because FFTW calculates an unnormalized DFT +constexpr auto scaleFactor = -float(M_GRAV) / (NCELL * 4); + +static_assert(sizeof(std::complex) == sizeof(fftwf_complex)); +struct FftwArrayDeleter { void operator ()(float ptr[]) const { fftwf_free(ptr); } }; +struct FftwComplexArrayDeleter { void operator ()(std::complex ptr[]) const { fftwf_free(ptr); } }; +struct FftwPlanDeleter { void operator ()(fftwf_plan ptr ) const { fftwf_destroy_plan(ptr); } }; +using FftwArrayPtr = std::unique_ptr; +using FftwComplexArrayPtr = std::unique_ptr [], FftwComplexArrayDeleter>; +using FftwPlanPtr = std::unique_ptr::type, FftwPlanDeleter >; +FftwArrayPtr FftwArray(size_t size) +{ + return FftwArrayPtr(reinterpret_cast(fftwf_malloc(size * sizeof(float)))); +} +FftwComplexArrayPtr FftwComplexArray(size_t size) +{ + return FftwComplexArrayPtr(reinterpret_cast *>(fftwf_malloc(size * sizeof(std::complex)))); +} + struct GravityImpl : public Gravity { bool grav_fft_status = false; - float *th_ptgravx = nullptr; - float *th_ptgravy = nullptr; - float *th_gravmapbig = nullptr; - float *th_gravxbig = nullptr; - float *th_gravybig = nullptr; - - fftwf_complex *th_ptgravxt, *th_ptgravyt, *th_gravmapbigt, *th_gravxbigt, *th_gravybigt; - fftwf_plan plan_gravmap, plan_gravx_inverse, plan_gravy_inverse; + FftwArrayPtr th_gravmapbig , th_gravxbig , th_gravybig ; + FftwComplexArrayPtr th_ptgravxt, th_ptgravyt, th_gravmapbigt, th_gravxbigt, th_gravybigt; + FftwPlanPtr plan_gravmap, plan_gravx_inverse, plan_gravy_inverse; void grav_fft_init(); void grav_fft_cleanup(); @@ -34,34 +52,28 @@ GravityImpl::~GravityImpl() void GravityImpl::grav_fft_init() { - int xblock2 = XCELLS*2; - int yblock2 = YCELLS*2; - int fft_tsize = (xblock2/2+1)*yblock2; - float distance, scaleFactor; - fftwf_plan plan_ptgravx, plan_ptgravy; if (grav_fft_status) return; + FftwPlanPtr plan_ptgravx, plan_ptgravy; //use fftw malloc function to ensure arrays are aligned, to get better performance - th_ptgravx = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_ptgravy = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_ptgravxt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); - th_ptgravyt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); - th_gravmapbig = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_gravmapbigt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); - th_gravxbig = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_gravybig = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_gravxbigt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); - th_gravybigt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); + FftwArrayPtr th_ptgravx = FftwArray(xblock2 * yblock2); + FftwArrayPtr th_ptgravy = FftwArray(xblock2 * yblock2); + th_ptgravxt = FftwComplexArray(fft_tsize); + th_ptgravyt = FftwComplexArray(fft_tsize); + th_gravmapbig = FftwArray(xblock2 * yblock2); + th_gravmapbigt = FftwComplexArray(fft_tsize); + th_gravxbig = FftwArray(xblock2 * yblock2); + th_gravybig = FftwArray(xblock2 * yblock2); + th_gravxbigt = FftwComplexArray(fft_tsize); + th_gravybigt = FftwComplexArray(fft_tsize); //select best algorithm, could use FFTW_PATIENT or FFTW_EXHAUSTIVE but that increases the time taken to plan, and I don't see much increase in execution speed - plan_ptgravx = fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_ptgravx, th_ptgravxt, FFTW_MEASURE); - plan_ptgravy = fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_ptgravy, th_ptgravyt, FFTW_MEASURE); - plan_gravmap = fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_gravmapbig, th_gravmapbigt, FFTW_MEASURE); - plan_gravx_inverse = fftwf_plan_dft_c2r_2d(yblock2, xblock2, th_gravxbigt, th_gravxbig, FFTW_MEASURE); - plan_gravy_inverse = fftwf_plan_dft_c2r_2d(yblock2, xblock2, th_gravybigt, th_gravybig, FFTW_MEASURE); + plan_ptgravx = FftwPlanPtr(fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_ptgravx.get(), reinterpret_cast(th_ptgravxt.get()), FFTW_MEASURE)); + plan_ptgravy = FftwPlanPtr(fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_ptgravy.get(), reinterpret_cast(th_ptgravyt.get()), FFTW_MEASURE)); + plan_gravmap = FftwPlanPtr(fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_gravmapbig.get(), reinterpret_cast(th_gravmapbigt.get()), FFTW_MEASURE)); + plan_gravx_inverse = FftwPlanPtr(fftwf_plan_dft_c2r_2d(yblock2, xblock2, reinterpret_cast(th_gravxbigt.get()), th_gravxbig.get(), FFTW_MEASURE)); + plan_gravy_inverse = FftwPlanPtr(fftwf_plan_dft_c2r_2d(yblock2, xblock2, reinterpret_cast(th_gravybigt.get()), th_gravybig.get(), FFTW_MEASURE)); - //NCELL*4 is size of data array, scaling needed because FFTW calculates an unnormalized DFT - scaleFactor = -float(M_GRAV)/(NCELL*4); //calculate velocity map caused by a point mass for (int y = 0; y < yblock2; y++) { @@ -69,7 +81,7 @@ void GravityImpl::grav_fft_init() { if (x == XCELLS && y == YCELLS) continue; - distance = hypotf(float(x-XCELLS), float(y-YCELLS)); + auto distance = hypotf(float(x-XCELLS), float(y-YCELLS)); th_ptgravx[y * xblock2 + x] = scaleFactor * (x - XCELLS) / powf(distance, 3); th_ptgravy[y * xblock2 + x] = scaleFactor * (y - YCELLS) / powf(distance, 3); } @@ -78,15 +90,11 @@ void GravityImpl::grav_fft_init() th_ptgravy[yblock2 * xblock2 / 2 + xblock2 / 2] = 0.0f; //transform point mass velocity maps - fftwf_execute(plan_ptgravx); - fftwf_execute(plan_ptgravy); - fftwf_destroy_plan(plan_ptgravx); - fftwf_destroy_plan(plan_ptgravy); - fftwf_free(th_ptgravx); - fftwf_free(th_ptgravy); + fftwf_execute(plan_ptgravx.get()); + fftwf_execute(plan_ptgravy.get()); //clear padded gravmap - memset(th_gravmapbig, 0, xblock2 * yblock2 * sizeof(float)); + memset(th_gravmapbig.get(), 0, xblock2 * yblock2 * sizeof(float)); grav_fft_status = true; } @@ -94,17 +102,6 @@ void GravityImpl::grav_fft_init() void GravityImpl::grav_fft_cleanup() { if (!grav_fft_status) return; - fftwf_free(th_ptgravxt); - fftwf_free(th_ptgravyt); - fftwf_free(th_gravmapbig); - fftwf_free(th_gravmapbigt); - fftwf_free(th_gravxbig); - fftwf_free(th_gravybig); - fftwf_free(th_gravxbigt); - fftwf_free(th_gravybigt); - fftwf_destroy_plan(plan_gravmap); - fftwf_destroy_plan(plan_gravx_inverse); - fftwf_destroy_plan(plan_gravy_inverse); grav_fft_status = false; } @@ -121,26 +118,23 @@ void Gravity::update_grav() if (!fftGravity->grav_fft_status) fftGravity->grav_fft_init(); - auto *th_gravmapbig = fftGravity->th_gravmapbig; - auto *th_gravxbig = fftGravity->th_gravxbig; - auto *th_gravybig = fftGravity->th_gravybig; - auto *th_ptgravxt = fftGravity->th_ptgravxt; - auto *th_ptgravyt = fftGravity->th_ptgravyt; - auto *th_gravmapbigt = fftGravity->th_gravmapbigt; - auto *th_gravxbigt = fftGravity->th_gravxbigt; - auto *th_gravybigt = fftGravity->th_gravybigt; + auto *th_gravmapbig = fftGravity->th_gravmapbig.get(); + auto *th_gravxbig = fftGravity->th_gravxbig.get(); + auto *th_gravybig = fftGravity->th_gravybig.get(); + auto *th_ptgravxt = fftGravity->th_ptgravxt.get(); + auto *th_ptgravyt = fftGravity->th_ptgravyt.get(); + auto *th_gravmapbigt = fftGravity->th_gravmapbigt.get(); + auto *th_gravxbigt = fftGravity->th_gravxbigt.get(); + auto *th_gravybigt = fftGravity->th_gravybigt.get(); auto &plan_gravmap = fftGravity->plan_gravmap; auto &plan_gravx_inverse = fftGravity->plan_gravx_inverse; auto &plan_gravy_inverse = fftGravity->plan_gravy_inverse; - int xblock2 = XCELLS*2, yblock2 = YCELLS*2; - int fft_tsize = (xblock2/2+1)*yblock2; - float mr, mc, pr, pc, gr, gc; - if (memcmp(th_ogravmap, th_gravmap, sizeof(float)*NCELL) != 0) + if (memcmp(&th_ogravmap[0], &th_gravmap[0], sizeof(float) * NCELL) != 0) { th_gravchanged = 1; - membwand(th_gravmap, gravmask, NCELL*sizeof(float), NCELL*sizeof(unsigned)); + membwand(&th_gravmap[0], &gravmask[0], NCELL * sizeof(float), NCELL * sizeof(uint32_t)); //copy gravmap into padded gravmap array for (int y = 0; y < YCELLS; y++) { @@ -150,28 +144,16 @@ void Gravity::update_grav() } } //transform gravmap - fftwf_execute(plan_gravmap); + fftwf_execute(plan_gravmap.get()); //do convolution (multiply the complex numbers) for (int i = 0; i < fft_tsize; i++) { - mr = th_gravmapbigt[i][0]; - mc = th_gravmapbigt[i][1]; - pr = th_ptgravxt[i][0]; - pc = th_ptgravxt[i][1]; - gr = mr*pr-mc*pc; - gc = mr*pc+mc*pr; - th_gravxbigt[i][0] = gr; - th_gravxbigt[i][1] = gc; - pr = th_ptgravyt[i][0]; - pc = th_ptgravyt[i][1]; - gr = mr*pr-mc*pc; - gc = mr*pc+mc*pr; - th_gravybigt[i][0] = gr; - th_gravybigt[i][1] = gc; + th_gravxbigt[i] = th_gravmapbigt[i] * th_ptgravxt[i]; + th_gravybigt[i] = th_gravmapbigt[i] * th_ptgravyt[i]; } //inverse transform, and copy from padded arrays into normal velocity maps - fftwf_execute(plan_gravx_inverse); - fftwf_execute(plan_gravy_inverse); + fftwf_execute(plan_gravx_inverse.get()); + fftwf_execute(plan_gravy_inverse.get()); for (int y = 0; y < YCELLS; y++) { for (int x = 0; x < XCELLS; x++) diff --git a/src/simulation/gravity/Gravity.h b/src/simulation/gravity/Gravity.h index 7e35d232e..bec94ab68 100644 --- a/src/simulation/gravity/Gravity.h +++ b/src/simulation/gravity/Gravity.h @@ -6,6 +6,8 @@ #include #include #include +#include +#include class Simulation; @@ -15,11 +17,11 @@ protected: bool enabled = false; // Maps to be processed by the gravity thread - float *th_ogravmap = nullptr; - float *th_gravmap = nullptr; - float *th_gravx = nullptr; - float *th_gravy = nullptr; - float *th_gravp = nullptr; + std::vector th_ogravmap; + std::vector th_gravmap; + std::vector th_gravx; + std::vector th_gravy; + std::vector th_gravp; int th_gravchanged = 0; @@ -53,11 +55,12 @@ public: ~Gravity(); //Maps to be used by the main thread - float *gravmap = nullptr; - float *gravp = nullptr; - float *gravy = nullptr; - float *gravx = nullptr; - unsigned *gravmask = nullptr; + std::vector gravmap; + std::vector gravp; + std::vector gravy; + std::vector gravx; + std::vector gravmask; + static_assert(sizeof(float) == sizeof(uint32_t)); unsigned char (*bmap)[XCELLS]; diff --git a/src/simulation/gravity/Null.cpp b/src/simulation/gravity/Null.cpp index 874d45eb4..02571fff6 100644 --- a/src/simulation/gravity/Null.cpp +++ b/src/simulation/gravity/Null.cpp @@ -6,9 +6,9 @@ void Gravity::get_result() { - memcpy(gravy, th_gravy, NCELL*sizeof(float)); - memcpy(gravx, th_gravx, NCELL*sizeof(float)); - memcpy(gravp, th_gravp, NCELL*sizeof(float)); + memcpy(&gravy[0], &th_gravy[0], NCELL * sizeof(float)); + memcpy(&gravx[0], &th_gravx[0], NCELL * sizeof(float)); + memcpy(&gravp[0], &th_gravp[0], NCELL * sizeof(float)); } void Gravity::update_grav(void)