Rebased splice potential

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

genetIC/src/ic.hpp

@@ -63,7 +63,7 @@ class ICGenerator {
   //! Vector of output fields (defaults to just a single field)
   std::vector<std::shared_ptr<fields::OutputField<GridDataType>>> outputFields;
   bool useBaryonTransferFunction{false}; //!< True if gas particles should use a different transfer function
-  modifications::ModificationManager<GridDataType> modificationManager; //!< Handles applying modificaitons to the various fields.
+  modifications::ModificationManager<GridDataType> modificationManager; //!< Handles applying modifications to the various fields.
   std::unique_ptr<fields::RandomFieldGenerator<GridDataType>> randomFieldGenerator; //!< Generate white noise for the output fields
   std::unique_ptr<cosmology::PowerSpectrum<GridDataType>> 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<particle::mapper::ParticleMapper<GridDataType>> 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; level<multiLevelContext.getNumLevels(); ++level) {
       auto &originalFieldThisLevel = outputFields[0]->getFieldForLevel(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) {
 

genetIC/src/main.cpp

@@ -137,16 +137,26 @@ void setup_parser(tools::ClassDispatch<ICType, void> &dispatch) {
 
   dispatch.add_class_route("reverse", static_cast<void (ICType::*)()>(&ICType::reverse));
   dispatch.add_class_route("reverse_small_k", static_cast<void (ICType::*)(FloatType)>(&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<void (ICType::*)(size_t)>(&ICType::dumpGrid));
   dispatch.add_class_route("dump_vx", static_cast<void (ICType::*)(size_t)>(&ICType::dumpVelocityX));
   dispatch.add_class_route("dump_grid_for_field", static_cast<void (ICType::*)(size_t, particle::species)>(&ICType::dumpGrid));
   dispatch.add_class_route("dump_grid_fourier", static_cast<void (ICType::*)(size_t)>(&ICType::dumpGridFourier));
-  dispatch.add_class_route("dump_grid_fourier_for_field",
-                           static_cast<void (ICType::*)(size_t, particle::species)>(&ICType::dumpGridFourier));
+  dispatch.add_class_route("dump_grid_fourier_for_field", static_cast<void (ICType::*)(size_t, particle::species)>(&ICType::dumpGridFourier));
   dispatch.add_class_route("dump_ps", static_cast<void (ICType::*)(size_t)>(&ICType::dumpPS));
   dispatch.add_class_route("dump_ps_field", static_cast<void (ICType::*)(size_t, particle::species)>(&ICType::dumpPS));
   dispatch.add_class_route("dump_tipsy", static_cast<void (ICType::*)(std::string)>(&ICType::saveTipsyArray));

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<DataType> max=0;
+      size_t N = data.size();
+
+      for(size_t i=0; i<N; i++) {
+        max = std::max(max, std::abs(data[i]));
+      }
+      return max;
+    }
+
     //! Single grid inner product. Only implemented for real space.
     auto innerProduct(const Field<DataType, CoordinateType> & other) const {
       assert(!other.isFourier());

genetIC/src/simulation/modifications/splice.hpp

@@ -4,6 +4,7 @@
 #include <complex>
 #include <src/tools/data_types/complex.hpp>
 #include <src/tools/numerics/cg.hpp>
+#include <src/tools/numerics/minres.hpp>
 
 namespace modifications {
   template<typename T>
@@ -37,7 +38,14 @@ namespace modifications {
   template<typename DataType, typename T=tools::datatypes::strip_complex<DataType>>
   fields::Field<DataType,T> spliceOneLevel(fields::Field<DataType,T> & a,
                                            fields::Field<DataType,T> & b,
-                                           const fields::Field<DataType,T> & cov) {
+                                           const fields::Field<DataType,T> & 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<DataType,T> preconditioner(cov);
+
       preconditioner.setFourierCoefficient(0, 0, 0, 1);
+      if (k_factor != 0) {
+        auto divide_by_k = [k_factor](std::complex<DataType> val, DataType kx, DataType ky, DataType kz){
+          DataType k2 = kx * kx + ky * ky + kz * kz;
+          if (k2 == 0 && k_factor < 0) {
+            return std::complex<DataType>(0, 0);
+          } else {
+            return val * DataType(pow(k2, k_factor));
+          }
+        };
+        preconditioner.forEachFourierCell(divide_by_k);
+      }
 
       fields::Field<DataType,T> delta_diff = b-a;
       delta_diff.applyTransferFunction(preconditioner, 0.5);
@@ -91,8 +111,14 @@ namespace modifications {
         return v;
       };
 
-
-      fields::Field<DataType,T> alpha = tools::numerics::conjugateGradient<DataType>(X,z);
+      fields::Field<DataType,T> alpha(z);
+      if (minimization_method == "CG") {
+        alpha = tools::numerics::conjugateGradient<DataType>(X, z, rtol, atol);
+      } else if (minimization_method == "MINRES") {
+        alpha = tools::numerics::minres<DataType>(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);

genetIC/src/tools/numerics/cg.hpp

@@ -5,6 +5,10 @@
 #include <src/simulation/field/field.hpp>
 #include <src/tools/data_types/complex.hpp>
 #include <src/tools/logging.hpp>
+#include <ctime>
+#include <iostream>
+#include <fstream>
+#include <sys/stat.h>
 
 namespace tools {
   namespace numerics {
@@ -14,18 +18,22 @@ namespace tools {
     fields::Field<T> conjugateGradient(std::function<fields::Field<T>(const fields::Field<T> &)> Q,
                                        const fields::Field<T> &b,
                                        double rtol = 1e-6,
-                                       double atol = 1e-12) {
+                                       double atol = 1e-12)
+    {
       fields::Field<T> residual(b);
       fields::Field<T> direction = -residual;
       fields::Field<T> x = fields::Field<T>(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;
-
     }
   }
 }