diff --git a/genetIC/Makefile b/genetIC/Makefile index f5c2bcd9..aefbd102 100644 --- a/genetIC/Makefile +++ b/genetIC/Makefile @@ -8,7 +8,7 @@ CFLAGS ?= -Wall -g -lpthread -O3 -fopenmp -std=c++17 -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 @@ -20,6 +20,7 @@ CFLAGS ?= -Wall -g -lpthread -O3 -fopenmp -std=c++17 -fdiagnostics-color=auto -I # -DFILTER_ON_COARSE_GRID: When moving long wavelength information onto a zoom grid, filters on the base grid # and then interpolates into the zoom, rather than the default which is to interpolate and then filter. +# Comment out -DZELDOVICH_GRADIENT_FOURIER_SPACE for potential splicing CODEOPTIONS = -DDOUBLEPRECISION -DOUTPUT_IN_DOUBLEPRECISION -DCUBIC_INTERPOLATION -DFRACTIONAL_K_SPLIT=0.3 -DFRACTIONAL_FILTER_TEMPERATURE=0.1 -DFILTER_ON_COARSE_GRID -DZELDOVICH_GRADIENT_FOURIER_SPACE CPATH ?= /opt/local/include/ LPATH ?= /opt/local/lib @@ -68,6 +69,14 @@ ifeq ($(USE_CUFFT), 1) CFLAGS=-O3 -I`pwd` -I`pwd`/HighFive/include -Xcompiler="-O3 -fopenmp -std=c++14 -Wall" CXX=nvcc endif + +ifeq ($(HOST3), Ana) + CXX = /opt/homebrew/bin/c++-14 + 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) 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 73a90960..e35eaa03 100644 --- a/genetIC/src/ic.hpp +++ b/genetIC/src/ic.hpp @@ -63,7 +63,7 @@ class ICGenerator { //! Vector of output fields (defaults to just a single field) std::vector>> outputFields; bool useBaryonTransferFunction{false}; //!< True if gas particles should use a different transfer function - modifications::ModificationManager modificationManager; //!< Handles applying modificaitons to the various fields. + modifications::ModificationManager modificationManager; //!< Handles applying modifications to the various fields. std::unique_ptr> randomFieldGenerator; //!< Generate white noise for the output fields std::unique_ptr> spectrum; //!< Transfer function data @@ -160,6 +160,27 @@ 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 minimization method for splicing + T splicing_rel_tol = 1e-6; + T splicing_abs_tol = 1e-12; + + //! Set the standard minimization method for splicing (CG or MINRES) + std::string minimization_method = "CG"; + + //! (does not work yet) Allows to continue minimizing a spliced field from a previous run + bool restart = false; + + //! Using fourier parallel seeding for splicing + bool setSplicedSeedFourierParallel = false; + //! Using fourier series seeding for splicing + bool setSplicedSeedFourierSeries = false; + //! Using real seeding for splicing + std::string spliceSeedingType = "Real"; + + //! The brake time in hours, set to 0 by default. If set to a positive value, splicing will stop after the given time, + // if it has not yet met the minimization threshold. + double brake_time = 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. @@ -259,7 +280,7 @@ class ICGenerator { allowStrayParticles = true; } -//!\brief Enables the use of a baryon density field. + //!\brief Enables the use of a baryon density field. //! If flagged on, the code will extract the transfer function for baryons from CAMB separately from //! the dark matter transfer. This causes twice as many fields to be stored and generated, but //! produces more accurate results for baryons than assuming they follow the same transfer function @@ -713,11 +734,6 @@ class ICGenerator { this->updateParticleMapper(); } - - - - - //! \brief Obtain power spectrum from a CAMB data file /*! * \param cambFieldPath - string of path to CAMB file @@ -1655,7 +1671,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 spliceWithFactor(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"); @@ -1673,19 +1689,87 @@ class ICGenerator { logging::entry() << "Constructing new random field for exterior of splice" << endl; newGenerator.seed(newSeed); + if(setSplicedSeedFourierParallel == true || setSplicedSeedFourierSeries == true) { + logging::entry() << "-!- Drawing splice seed in " << spliceSeedingType << " -!-" << endl; + newGenerator.setDrawInFourierSpace(true); + if(setSplicedSeedFourierParallel == 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_rel_tol, + splicing_abs_tol, + k_factor, + minimization_method, + restart, + brake_time + ); 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 setSpliceAccuracy(string type, double tolerance) { + if (type == "absolute") { + splicing_abs_tol = tolerance; + logging::entry() << "Setting splicing " << minimization_method << " absolute tolerance to " << tolerance << std::endl; + } else if (type == "relative") { + splicing_rel_tol = tolerance; + logging::entry() << "Setting splicing " << minimization_method << " relative tolerance to " << tolerance << std::endl; + } + } + + virtual void spliceSeedFourierSeries() { + setSplicedSeedFourierSeries = true; + spliceSeedingType = "fourier and series"; + } + + virtual void spliceSeedFourierParallel() { + setSplicedSeedFourierParallel = true; + spliceSeedingType = "fourier and parallel"; + } + + virtual void restartSplice() { + restart = true; + } + + virtual void setSpliceMinimization(string set_minMethod) { + minimization_method = set_minMethod; + } + + virtual void stopSplicingAfter(double time_to_brake) { + brake_time = time_to_brake; + } + + virtual void spliceDensity(size_t newSeed) { + spliceWithFactor(newSeed, 0); + } + + virtual void splicePotential(size_t newSeed) { + #ifdef ZELDOVICH_GRADIENT_FOURIER_SPACE + logging::entry(logging::level::warning) << "WARNING: You are splicing the potential field with ZELDOVICH_GRADIENT_FOURIER_SPACE." << std::endl; + logging::entry(logging::level::warning) << " Due to non-locality, the density field may not exactly match within the" << std::endl; + logging::entry(logging::level::warning) << " spliced region. Consider recompiling with ZELDOVICH_GRADIENT_FOURIER_SPACE" << std::endl; + logging::entry(logging::level::warning) << " disabled if precise density field matching is important to you." << std::endl; + #endif + spliceWithFactor(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 f5ef40fd..e28d3209 100644 --- a/genetIC/src/main.cpp +++ b/genetIC/src/main.cpp @@ -137,16 +137,26 @@ 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::setSpliceAccuracy); + dispatch.add_class_route("set_splice_accuracy", &ICType::setSpliceAccuracy); + dispatch.add_class_route("splice_seedfourier_parallel", &ICType::spliceSeedFourierParallel); + dispatch.add_class_route("splice_seedfourier_series", &ICType::spliceSeedFourierSeries); + dispatch.add_deprecated_class_route("set_minimization", "set_splice_minimization", &ICType::setSpliceMinimization); + dispatch.add_class_route("set_splice_minimization", &ICType::setSpliceMinimization); + dispatch.add_deprecated_class_route("splice", "splice_density", &ICType::spliceDensity); + dispatch.add_class_route("splice_density", &ICType::spliceDensity); + dispatch.add_class_route("splice_potential", &ICType::splicePotential); + + dispatch.add_deprecated_class_route("stop_after", "stop_splicing_after", &ICType::stopSplicingAfter); + dispatch.add_class_route("stop_splicing_after", &ICType::stopSplicingAfter); + // Write objects to files // dispatch.add_class_route("dump_grid", &ICType::dumpGrid); dispatch.add_class_route("dump_grid", static_cast(&ICType::dumpGrid)); dispatch.add_class_route("dump_vx", static_cast(&ICType::dumpVelocityX)); dispatch.add_class_route("dump_grid_for_field", static_cast(&ICType::dumpGrid)); dispatch.add_class_route("dump_grid_fourier", static_cast(&ICType::dumpGridFourier)); - dispatch.add_class_route("dump_grid_fourier_for_field", - static_cast(&ICType::dumpGridFourier)); + dispatch.add_class_route("dump_grid_fourier_for_field", static_cast(&ICType::dumpGridFourier)); dispatch.add_class_route("dump_ps", static_cast(&ICType::dumpPS)); dispatch.add_class_route("dump_ps_field", static_cast(&ICType::dumpPS)); dispatch.add_class_route("dump_tipsy", static_cast(&ICType::saveTipsyArray)); diff --git a/genetIC/src/simulation/field/field.hpp b/genetIC/src/simulation/field/field.hpp index 08f8659b..2b87a4b7 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..b2cc60b2 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 @@ -37,7 +38,14 @@ 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 double brake_time=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 @@ -55,7 +63,19 @@ namespace modifications { // We however set the fundamental of the power spectrum to a non-null value, // 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); @@ -91,8 +111,14 @@ namespace modifications { return v; }; - - fields::Field alpha = tools::numerics::conjugateGradient(X,z); + fields::Field alpha(z); + if (minimization_method == "CG") { + alpha = tools::numerics::conjugateGradient(X, z, rtol, atol); + } else if (minimization_method == "MINRES") { + alpha = tools::numerics::minres(X, z, rtol, atol, restart, brake_time); + } else { + throw std::runtime_error("Minimization method is invalid. Current implementations are CG and MINRES"); + } alpha.toFourier(); alpha.applyTransferFunction(preconditioner, 0.5); diff --git a/genetIC/src/tools/numerics/cg.hpp b/genetIC/src/tools/numerics/cg.hpp index bce1cf18..e555de58 100644 --- a/genetIC/src/tools/numerics/cg.hpp +++ b/genetIC/src/tools/numerics/cg.hpp @@ -5,6 +5,10 @@ #include #include #include +#include +#include +#include +#include namespace tools { namespace numerics { @@ -14,18 +18,22 @@ namespace tools { fields::Field conjugateGradient(std::function(const fields::Field &)> Q, const fields::Field &b, double rtol = 1e-6, - double atol = 1e-12) { + double atol = 1e-12) + { fields::Field residual(b); fields::Field direction = -residual; fields::Field x = fields::Field(b.getGrid(), false); - double scale = b.norm(); + // double scaleNorm = b.norm(); + double scaleMax = b.maximum(); - if(scale==0.0) { + if(scaleMax==0.0) { logging::entry(logging::warning) << "Conjugate gradient: result is zero!" << std::endl; return x; } + logging::entry() << "Splicing will stop when the maximum drops below " << rtol * scaleMax << std::endl; + size_t dimension = b.getGrid().size3; size_t i; @@ -40,11 +48,13 @@ namespace tools { residual = Q(x); residual -= b; - auto norm = residual.norm(); - if (norm < rtol * scale || norm < atol) + // auto norm = residual.norm(); + auto norm = residual.maximum(); + + if (norm < rtol * scaleMax || norm < atol) break; - logging::entry() << "Conjugate gradient iteration " << i << " residual=" << norm << std::endl; + logging::entry() << "CG iteration " << i << " maximum=" << norm << std::endl; // update direction for next cycle; must be Q-orthogonal to all previous updates double beta = residual.innerProduct(Q_direction) / direction.innerProduct(Q_direction); @@ -55,7 +65,6 @@ namespace tools { logging::entry() << "Conjugate gradient ended after " << i << " iterations" << std::endl; return x; - } } } diff --git a/genetIC/src/tools/numerics/minres.hpp b/genetIC/src/tools/numerics/minres.hpp new file mode 100644 index 00000000..8904aff0 --- /dev/null +++ b/genetIC/src/tools/numerics/minres.hpp @@ -0,0 +1,210 @@ +#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 brake_time = 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; + brake_time = 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; + + // Set variables from last restart (to continue minimization on a spliced field from a previous run) + // -! not working yet !- + /* + 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 (brake_time > 0) { + + const std::time_t brakeDuration = brake_time * 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 " << brake_time << " hours, or when the maximum drops below " << rtol * scaleMax << std::endl; + + // Start iteration + for(; iter