Splice potential

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/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

@@ -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<particle::mapper::ParticleMapper<GridDataType>> 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; 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_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) {
 

genetIC/src/main.cpp

@@ -135,7 +135,16 @@ 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::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);

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,15 @@ 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 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 +65,17 @@ namespace modifications {
       // 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);
@@ -92,29 +112,57 @@ namespace modifications {
       };
 
 
-      fields::Field<DataType,T> alpha = tools::numerics::conjugateGradient<DataType>(X,z);
-
-      alpha.toFourier();
-      alpha.applyTransferFunction(preconditioner, 0.5);
-      alpha.toReal();
-
-      fields::Field<DataType,T> 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<DataType,T> alpha = tools::numerics::conjugateGradient<DataType>(X, z, rtol, atol);
+        alpha.toFourier();
+        alpha.applyTransferFunction(preconditioner, 0.5);
+        alpha.toReal();
+
+        fields::Field<DataType,T> 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<DataType,T> alpha = tools::numerics::minres<DataType>(X, z, rtol, atol, restart, stop, brakeTime);
+        alpha.toFourier();
+        alpha.applyTransferFunction(preconditioner, 0.5);
+        alpha.toReal();
+
+        fields::Field<DataType,T> 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");
+      }
   }
 }