From d8f9ebe98f9716ad34959f50fc14e204d000eb86 Mon Sep 17 00:00:00 2001 From: Anatole Storck Date: Fri, 10 Nov 2023 18:48:32 +0100 Subject: [PATCH 1/4] Adding potential splicing capabilities. Along with a new minimization method (MINRES) which is much faster than the current CG method. Using MINRES is a must for potential splicing of large fields (scaling is much better). --- genetIC/Makefile | 10 +- genetIC/src/ic.hpp | 89 +++++++- genetIC/src/main.cpp | 11 +- genetIC/src/simulation/field/field.hpp | 13 ++ .../src/simulation/modifications/splice.hpp | 95 +++++--- genetIC/src/tools/numerics/minres.hpp | 202 ++++++++++++++++++ 6 files changed, 390 insertions(+), 30 deletions(-) create mode 100644 genetIC/src/tools/numerics/minres.hpp diff --git a/genetIC/Makefile b/genetIC/Makefile index 7a0df0ac..64f5ab5d 100644 --- a/genetIC/Makefile +++ b/genetIC/Makefile @@ -8,7 +8,7 @@ CFLAGS ?= -Wall -g -lpthread -O3 -fopenmp -std=c++14 -fdiagnostics-color=auto -I # -DVELOCITY_MODIFICATION_GRADIENT_FOURIER_SPACE: As the name suggests, use "ik" as the gradient operator; otherwise # uses a fourth order stencil in real space. Enabling this is essential if you disable CUBIC_INTERPOLATION. # -DZELDOVICH_GRADIENT_FOURIER_SPACE: Same as VELOCITY_MODIFICATION_GRADIENT_FOURIER_SPACE, but for computing the -# Zeldovich displacement field. +# Zeldovich displacement field. Should be commented out if splicing the potential. # -DFRACTIONAL_K_SPLIT: This defines the characteristic k0 of the filter used to split information between zoom levels, # defined by k0 = k_nyquist * FRACTIONAL_K_SPLIT, where k_nyquist is the Nyquist frequency of the lower-resolution grid. # Higher values correlate the zoom grids more precisely with the unzoomed base grids, but at the cost of errors @@ -61,6 +61,14 @@ ifeq ($(USE_CUFFT), 1) CFLAGS=-O3 -I`pwd` -Xcompiler="-O3 -fopenmp -std=c++14 -Wall" CXX=nvcc endif + +ifeq ($(HOST3), Ana) + CXX = /opt/homebrew/bin/c++-13 + CPATH ?= /opt/local/include/ + LPATH ?= /opt/local/lib + CFLAGS = -O3 -fopenmp -std=c++14 -fdiagnostics-color=auto -I`pwd` -DOPENMP -DDOUBLEPRECISION +endif + ifeq ($(HOST3), pfe) CXX = /nasa/pkgsrc/2016Q2/gcc5/bin/g++ CPATH = /u/apontzen/genetIC/genetIC:/nasa/gsl/1.14/include/:/nasa/intel/Compiler/2016.2.181/compilers_and_libraries_2016.2.181/linux/mkl/include/fftw diff --git a/genetIC/src/ic.hpp b/genetIC/src/ic.hpp index 9e7ac0f6..01f5f4ba 100644 --- a/genetIC/src/ic.hpp +++ b/genetIC/src/ic.hpp @@ -160,6 +160,25 @@ class ICGenerator { int initial_number_steps = 10; //!< Number of steps always used by quadratic modifications. T precision = 0.001; //!< Target precision required of quadratic modifications. + //! Accuracy of the conjugate gradient method for the splicing method + T splicing_cg_rel_tol = 1e-6; + T splicing_cg_abs_tol = 1e-12; + + //! Set the standard minimization method for splicing + std::string minimization_method = "CG"; + + //! Restarting function for splicing + bool restart = false; + + //! Using fourier parallel seeding for splicing + bool setSplicedSeedFourierParallel = false; + + //! Wether to stop splicing after a certain amount of time + bool stop = false; + + //! Set the brake time to zero (unused) + double brakeTime = 0; + //! Mapper that keep track of particles in the mulit-level context. shared_ptr> pMapper = nullptr; //! Input mapper, used to relate particles in a different simulation to particles in this one. @@ -1635,7 +1654,7 @@ class ICGenerator { } //! Splicing: fixes the flagged region, while reinitialising the exterior from a new random field - virtual void splice(size_t newSeed) { + virtual void splice_with_factor(size_t newSeed, int k_factor=0) { initialiseRandomComponentIfUninitialised(); if(outputFields.size()>1) throw std::runtime_error("Splicing is not yet implemented for the case of multiple transfer functions"); @@ -1653,19 +1672,81 @@ class ICGenerator { logging::entry() << "Constructing new random field for exterior of splice" << endl; newGenerator.seed(newSeed); + if(setSplicedSeedFourierParallel == true) { + logging::entry() << "-!- Drawing splice seed in fourier and parallel -!-" << endl; + newGenerator.setDrawInFourierSpace(true); + newGenerator.setParallel(true); + newGenerator.setReverseRandomDrawOrder(false); + } newGenerator.draw(); - logging::entry() << "Finished constructing new random field. Beginning splice operation." << endl; + logging::entry() << "Finished constructing new random field" << endl; + logging::entry() << "Beginning splice operation using the " << minimization_method << " method." << endl; for(size_t level=0; levelgetFieldForLevel(level); auto &newFieldThisLevel = newField.getFieldForLevel(level); - auto splicedFieldThisLevel = modifications::spliceOneLevel(newFieldThisLevel, originalFieldThisLevel, - *multiLevelContext.getCovariance(level, particle::species::all)); + auto splicedFieldThisLevel = modifications::spliceOneLevel( + newFieldThisLevel, + originalFieldThisLevel, + *multiLevelContext.getCovariance(level, particle::species::all), + splicing_cg_rel_tol, + splicing_cg_abs_tol, + k_factor, + minimization_method, + restart, + stop, + brakeTime + ); splicedFieldThisLevel.toFourier(); originalFieldThisLevel = std::move(splicedFieldThisLevel); } } + //! Set the conjugate gradient precision for the splicing method + /*! + * @param type Precision can be relative or absolute + * @param tolerance Precision of the CG + */ + virtual void set_splice_accuracy(string type, double tolerance) { + if (type == "absolute") { + splicing_cg_abs_tol = tolerance; + logging::entry() << "Setting splicing CG absolute tolerance to " << tolerance << std::endl; + } else if (type == "relative") { + splicing_cg_rel_tol = tolerance; + logging::entry() << "Setting splicing CG relative tolerance to " << tolerance << std::endl; + } + } + + virtual void splice_seedfourier_parallel() { + setSplicedSeedFourierParallel = true; + } + + virtual void restart_splice() { + restart = true; + } + + virtual void set_minimization(string set_minMethod) { + minimization_method = set_minMethod; + } + + virtual void stop_after(double time_to_brake) { + stop = true; + brakeTime = time_to_brake; + } + + virtual void splice_density(size_t newSeed) { + splice_with_factor(newSeed, 0); + } + + virtual void splice_velocity(size_t newSeed) { + splice_with_factor(newSeed, -1); + } + + virtual void splice_potential(size_t newSeed) { + // minimization_method = "MINRES"; + splice_with_factor(newSeed, -2); + } + //! Reverses the sign of the low-k modes. virtual void reverseSmallK(T kmax) { diff --git a/genetIC/src/main.cpp b/genetIC/src/main.cpp index db85b51f..83d6cd5e 100644 --- a/genetIC/src/main.cpp +++ b/genetIC/src/main.cpp @@ -135,7 +135,16 @@ void setup_parser(tools::ClassDispatch &dispatch) { dispatch.add_class_route("reverse", static_cast(&ICType::reverse)); dispatch.add_class_route("reverse_small_k", static_cast(&ICType::reverseSmallK)); - dispatch.add_class_route("splice", &ICType::splice); + dispatch.add_deprecated_class_route("splice_accuracy", "set_splice_accuracy", &ICType::set_splice_accuracy); + dispatch.add_class_route("set_splice_accuracy", &ICType::set_splice_accuracy); + dispatch.add_class_route("splice_seedfourier_parallel", &ICType::splice_seedfourier_parallel); + dispatch.add_class_route("restart_splice", &ICType::restart_splice); + dispatch.add_class_route("set_minimization", &ICType::set_minimization); + dispatch.add_deprecated_class_route("splice", "splice_density", &ICType::splice_density); + dispatch.add_class_route("splice_density", &ICType::splice_density); + dispatch.add_class_route("splice_potential", &ICType::splice_potential); + + dispatch.add_class_route("stop_after", &ICType::stop_after); // Write objects to files // dispatch.add_class_route("dump_grid", &ICType::dumpGrid); diff --git a/genetIC/src/simulation/field/field.hpp b/genetIC/src/simulation/field/field.hpp index 08f8659b..8c0d1135 100644 --- a/genetIC/src/simulation/field/field.hpp +++ b/genetIC/src/simulation/field/field.hpp @@ -361,6 +361,19 @@ namespace fields { } } + //! Single grid maximum. Only implemented for real space + auto Maximum() const { + assert(!isFourier()); + + tools::datatypes::strip_complex max=0; + size_t N = data.size(); + + for(size_t i=0; i & other) const { assert(!other.isFourier()); diff --git a/genetIC/src/simulation/modifications/splice.hpp b/genetIC/src/simulation/modifications/splice.hpp index 2cda17f4..994c0fb2 100644 --- a/genetIC/src/simulation/modifications/splice.hpp +++ b/genetIC/src/simulation/modifications/splice.hpp @@ -37,7 +37,15 @@ namespace modifications { template> fields::Field spliceOneLevel(fields::Field & a, fields::Field & b, - const fields::Field & cov) { + const fields::Field & cov, + const double rtol, + const double atol, + const int k_factor=0, + const std::string minimization_method="CG", + const bool restart=false, + const bool stop=false, + const double brakeTime=0 + ) { // To understand the implementation below, first read Appendix A of Cadiou et al (2021), // and/or look at the 1D toy implementation (in tools/toy_implementation/gene_splicing.ipynb) which @@ -56,6 +64,17 @@ namespace modifications { // otherwise, the mean value in the spliced region is unconstrained. fields::Field preconditioner(cov); preconditioner.setFourierCoefficient(0, 0, 0, 1); + if (k_factor != 0) { + auto divide_by_k = [k_factor](std::complex val, DataType kx, DataType ky, DataType kz){ + DataType k2 = kx * kx + ky * ky + kz * kz; + if (k2 == 0 && k_factor < 0) { + return std::complex(0, 0); + } else { + return val * DataType(pow(k2, k_factor)); + } + }; + preconditioner.forEachFourierCell(divide_by_k); + } fields::Field delta_diff = b-a; delta_diff.applyTransferFunction(preconditioner, 0.5); @@ -92,29 +111,57 @@ namespace modifications { }; - fields::Field alpha = tools::numerics::conjugateGradient(X,z); - - alpha.toFourier(); - alpha.applyTransferFunction(preconditioner, 0.5); - alpha.toReal(); - - fields::Field bInDeltaBasis(b); - bInDeltaBasis.toFourier(); - bInDeltaBasis.applyTransferFunction(preconditioner, 0.5); - bInDeltaBasis.toReal(); - - alpha*=maskCompl; - alpha+=bInDeltaBasis; - - delta_diff*=mask; - alpha-=delta_diff; - - assert(!alpha.isFourier()); - alpha.toFourier(); - alpha.applyTransferFunction(preconditioner, -0.5); - alpha.toReal(); - - return alpha; + if (minimization_method == "CG") { + fields::Field alpha = tools::numerics::conjugateGradient(X, z, rtol, atol); + alpha.toFourier(); + alpha.applyTransferFunction(preconditioner, 0.5); + alpha.toReal(); + + fields::Field bInDeltaBasis(b); + bInDeltaBasis.toFourier(); + bInDeltaBasis.applyTransferFunction(preconditioner, 0.5); + bInDeltaBasis.toReal(); + + alpha*=maskCompl; + alpha+=bInDeltaBasis; + + delta_diff*=mask; + alpha-=delta_diff; + + assert(!alpha.isFourier()); + alpha.toFourier(); + alpha.applyTransferFunction(preconditioner, -0.5); + alpha.toReal(); + + return alpha; + } + else if (minimization_method == "MINRES") { + fields::Field alpha = tools::numerics::minres(X, z, rtol, atol, restart, stop, brakeTime); + alpha.toFourier(); + alpha.applyTransferFunction(preconditioner, 0.5); + alpha.toReal(); + + fields::Field bInDeltaBasis(b); + bInDeltaBasis.toFourier(); + bInDeltaBasis.applyTransferFunction(preconditioner, 0.5); + bInDeltaBasis.toReal(); + + alpha*=maskCompl; + alpha+=bInDeltaBasis; + + delta_diff*=mask; + alpha-=delta_diff; + + assert(!alpha.isFourier()); + alpha.toFourier(); + alpha.applyTransferFunction(preconditioner, -0.5); + alpha.toReal(); + + return alpha; + } + else { + throw std::runtime_error("Minimization method is invalid. Current implementations are CG and MINRES"); + } } } diff --git a/genetIC/src/tools/numerics/minres.hpp b/genetIC/src/tools/numerics/minres.hpp new file mode 100644 index 00000000..19d67995 --- /dev/null +++ b/genetIC/src/tools/numerics/minres.hpp @@ -0,0 +1,202 @@ +#ifndef IC_MINRES_HPP +#define IC_MINRES_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace tools { + namespace numerics { + + /* Adapted from https://stanford.edu/group/SOL/reports/SOL-2011-2R.pdf */ + template + fields::Field minres(std::function(fields::Field &)> A, + fields::Field &b, + double rtol = 1e-6, + double atol = 1e-12, + bool restart = false, + bool stop = false, + double brakeTime = 0) + { + fields::Field x(b); + x *= 0; + fields::Field r(b); + fields::Field s(b); + s = A(r); + + fields::Field p(r); + fields::Field q(s); + + // auto toltest = 1e-8; // Default tolerance (meant to even include 1024^3 potential splicing) + + double scaleNorm = sqrt(b.innerProduct(b)); + double scaleMax = b.Maximum(); + double rho = r.innerProduct(s); + + size_t dimension = b.getGrid().size3; + + /* + if (rtol != 2e-4) { + + // Lowering tolerance for lower resolution fields + if (dimension == 512*512*512) + toltest = 1e-7; + if (dimension == 256*256*256) { + toltest = 1e-6; + brakeTime = 24; // Perhaps a time dependence instead of trying to achieve the + // exact tolerance is better + } + } + */ + + size_t max_iterations = dimension * 10; + double old_norm = 0; + + size_t iter = 0; + //setting variables from last restart + /* + if (restart == true) { + logging::entry() << "Loading restart variables" << std::endl; + + std::string variables = ".variables"; + std::string directory = output_path + variables; + + x.loadGridData(directory + "/x"); + r.loadGridData(directory + "/r"); + q.loadGridData(directory + "/q"); + p.loadGridData(directory + "/p"); + + std::ifstream file; + + file.open (directory + "/rho.txt"); + file >> rho; + file.close(); + + file.open (directory + "/iter.txt"); + file >> iter; + file.close(); + + } + */ + + if (stop == true){ + + const std::time_t brakeDuration = brakeTime * 3600; // Calculate the duration in seconds for the brake time + const std::time_t startTime = std::time(nullptr); // Get the current time at splicing start + + logging::entry() << "Splicing will stop after " << brakeTime << " hours, or when the maximum drops below " << rtol * scaleMax << std::endl; + + // Start iteration + for(; iter Date: Fri, 10 Nov 2023 18:51:17 +0100 Subject: [PATCH 2/4] need to include minres.hpp to splice.hpp --- genetIC/src/simulation/modifications/splice.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/genetIC/src/simulation/modifications/splice.hpp b/genetIC/src/simulation/modifications/splice.hpp index 994c0fb2..62c2fb86 100644 --- a/genetIC/src/simulation/modifications/splice.hpp +++ b/genetIC/src/simulation/modifications/splice.hpp @@ -4,6 +4,7 @@ #include #include #include +#include namespace modifications { template From 53ecb7e619e264d60c31d31bf5a312d8797e59c0 Mon Sep 17 00:00:00 2001 From: Anatole Storck Date: Fri, 10 Nov 2023 18:52:15 +0100 Subject: [PATCH 3/4] removed typo --- genetIC/src/tools/numerics/minres.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/genetIC/src/tools/numerics/minres.hpp b/genetIC/src/tools/numerics/minres.hpp index 19d67995..cf5d58bd 100644 --- a/genetIC/src/tools/numerics/minres.hpp +++ b/genetIC/src/tools/numerics/minres.hpp @@ -199,4 +199,4 @@ namespace tools { } } -+#endif \ No newline at end of file +#endif \ No newline at end of file From 64291fa54c17055c3f96102804a506a7c1d4fed8 Mon Sep 17 00:00:00 2001 From: Anatole Storck Date: Fri, 10 Nov 2023 18:57:32 +0100 Subject: [PATCH 4/4] updated paths for local machine --- genetIC/Makefile | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/genetIC/Makefile b/genetIC/Makefile index 64f5ab5d..f98e164c 100644 --- a/genetIC/Makefile +++ b/genetIC/Makefile @@ -63,10 +63,10 @@ ifeq ($(USE_CUFFT), 1) endif ifeq ($(HOST3), Ana) - CXX = /opt/homebrew/bin/c++-13 - CPATH ?= /opt/local/include/ - LPATH ?= /opt/local/lib - CFLAGS = -O3 -fopenmp -std=c++14 -fdiagnostics-color=auto -I`pwd` -DOPENMP -DDOUBLEPRECISION + CXX = /opt/homebrew/bin/c++-13 + CPATH = /opt/homebrew/include/ + LPATH = /opt/homebrew/lib + CFLAGS = -O3 -fopenmp -std=c++14 -fdiagnostics-color=auto -I`pwd` -DOPENMP -DDOUBLEPRECISION endif ifeq ($(HOST3), pfe)