[WIP] Splice multilevel

CMakeLists.txt

@@ -60,6 +60,7 @@ set(SOURCE_FILES
         genetIC/src/simulation/modifications/linearmodification.hpp
         genetIC/src/simulation/modifications/quadraticmodification.hpp
         genetIC/src/simulation/multilevelgrid/mask.hpp
+        genetIC/src/tools/numerics/cg.hpp
         genetIC/src/tools/memmap.hpp
         genetIC/src/tools/numerics/tricubic.hpp genetIC/src/tools/logging.hpp genetIC/src/tools/logging.cpp genetIC/src/simulation/modifications/splice.hpp genetIC/src/tools/lru_cache.hpp)
 
@@ -116,4 +117,9 @@ add_definitions(-DCUBIC_INTERPOLATION
         -DZELDOVICH_GRADIENT_FOURIER_SPACE
         )
 add_compile_options(-Wextra)
-add_executable(genetIC ${SOURCE_FILES})
+add_executable(genetIC ${SOURCE_FILES})
+find_package(OpenMP)
+if(OpenMP_CXX_FOUND)
+    target_link_libraries(genetIC PUBLIC OpenMP::OpenMP_CXX)
+    add_definitions(-DOPENMP)
+endif()

genetIC/src/ic.hpp

@@ -174,13 +174,16 @@ class ICGenerator {
 
   tools::ClassDispatch<ICGenerator<GridDataType>, void> &interpreter; //!< Parser for parameter files
 
+  bool shouldSplice; //!< True if the code should perform the splicing operation.
+  size_t spliceSeed; //!< Seed for the random number generator used in the splicing operation.
 
 public:
   //! \brief Main constructor for this class.
   //! Requires a single interpreter to initialise, from which it will receive input commands.
   ICGenerator(tools::ClassDispatch<ICGenerator<GridDataType>, void> &interpreter) :
     modificationManager(multiLevelContext, cosmology, nullptr),
-    interpreter(interpreter) {
+    interpreter(interpreter),
+    shouldSplice(false) {
 
     // By default, we assume there is only one field - at first this will contain white noise:
     outputFields.push_back(
@@ -849,7 +852,6 @@ class ICGenerator {
   */
   template<typename TField>
   void dumpGridData(size_t level, const TField &data, std::string prefix = "grid") {
-
     initialiseRandomComponentIfUninitialised();
 
     auto levelGrid = data.getGrid();
@@ -903,6 +905,7 @@ class ICGenerator {
   * \param species - the field to dump
   */
   virtual void dumpGrid(size_t level, particle::species species ) {
+    checkLevelExists(level, species);
     auto & field = this->getOutputFieldForSpecies(species);
     field.toReal();
     dumpGridData(level, field.getFieldForLevel(level));
@@ -917,17 +920,20 @@ class ICGenerator {
   }
 
   virtual void dumpVelocityX(size_t level) {
+    checkLevelExists(level, particle::species::dm);
     dumpGridData(level, *(this->pParticleGenerator[particle::species::dm]->getGeneratorForLevel(
       level).getGeneratedFields()[0]), "vx");
   }
 
   //! For backwards compatibility. Dumpts baryons to field at requested level to file named grid-level.
   virtual void dumpGridFourier(size_t level = 0) {
+    checkLevelExists(level, particle::species::dm);
     this->dumpGridFourier(level, particle::species::dm);
   }
 
   //! Output the grid in Fourier space.
   virtual void dumpGridFourier(size_t level, particle::species species) {
+    checkLevelExists(level, species);
     auto & field = this->getOutputFieldForSpecies(species);
     field.toFourier();
     fields::Field<complex<T>, T> fieldToWrite = tools::numerics::fourier::getComplexFourierField(field.getFieldForLevel(level));
@@ -1634,9 +1640,12 @@ class ICGenerator {
     outputFields[0]->reverse();
   }
 
-  //! Splicing: fixes the flagged region, while reinitialising the exterior from a new random field
+  //! Save that splicing should happen in the generation of the output
   virtual void splice(size_t newSeed) {
+    this->spliceSeed = newSeed;
+    this->shouldSplice = true;
     initialiseRandomComponentIfUninitialised();
+
     if(outputFields.size()>1)
       throw std::runtime_error("Splicing is not yet implemented for the case of multiple transfer functions");
 
@@ -1645,8 +1654,14 @@ class ICGenerator {
       throw std::runtime_error("It is too late in the IC generation process to perform splicing; try moving the splice command earlier");
     }
 
-    fields::OutputField<GridDataType> newField = fields::OutputField<GridDataType>(multiLevelContext, particle::species::whitenoise);
-    auto newGenerator = fields::RandomFieldGenerator<GridDataType>(newField);
+    spliceFields(newSeed);
+  }
+
+  //! Splicing: fixes the flagged region, while reinitialising the exterior from a new random field
+  virtual void spliceFields(size_t newSeed) {
+
+    fields::OutputField<GridDataType> a = fields::OutputField<GridDataType>(multiLevelContext, particle::species::whitenoise);
+    auto newGenerator = fields::RandomFieldGenerator<GridDataType>(a);
 
     if(multiLevelContext.getNumLevels()>1)
       logging::entry(logging::level::warning) << "Splicing operations on multiple levels work, but may have edge artefacts. Use experimentally/at your own risk!" << endl;
@@ -1655,15 +1670,20 @@ class ICGenerator {
     newGenerator.seed(newSeed);
     newGenerator.draw();
     logging::entry() << "Finished constructing new random field. Beginning splice operation." << endl;
+
+    auto & b = *outputFields[0];
 
-    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));
-      splicedFieldThisLevel.toFourier();
-      originalFieldThisLevel = std::move(splicedFieldThisLevel);
-    }
+    fields::OutputField<GridDataType> f = modifications::spliceMultiLevel<GridDataType>(b, a);
+
+    // Copy the result back into the output field
+    outputFields[0]->copyData(f);
+
+    // The field f is in delta space and we just copied it, so mark
+    // the outputFields accordingly. 
+    // outputFields[0]->setForceTransferType(particle::species::all);
+
+    // Mark the field as already combined
+    // multiLevelContext.setLevelsAreCombined();
   }
 
   //! Reverses the sign of the low-k modes.

genetIC/src/simulation/field/field.hpp

@@ -354,6 +354,9 @@ namespace fields {
     //! Add a multiple of the provided field to this one in-place
     void addScaled(const Field<DataType, CoordinateType> & other,
                    tools::datatypes::strip_complex<DataType> scale) {
+      assert(!other.isFourier());
+      assert(!isFourier());
+
       size_t N = data.size();
 #pragma omp parallel for
       for(size_t i=0; i<N; i++) {
@@ -367,7 +370,7 @@ namespace fields {
       assert(!isFourier());
 
       tools::datatypes::strip_complex<DataType> v=0;
-      size_t N = data.size();
+      size_t N = getGrid().size3;
 
 #pragma omp parallel for reduction(+:v)
       for(size_t i=0; i<N; i++) {
@@ -479,6 +482,35 @@ namespace fields {
 
     }
 
+    void deInterpolate(const Field<DataType, CoordinateType> &hires) {
+      assert (!isFourier());
+      assert (!hires.isFourier());
+      int pixel_size_ratio = tools::getRatioAndAssertInteger(getGrid().cellSize,
+                                                             hires.getGrid().cellSize);
+      int pixel_volume_ratio = pow(pixel_size_ratio,3);
+
+      size_t hiresGridSize = hires.getGrid().size;
+      size_t hiresParallelisationSeparation = 4 * pixel_size_ratio;
+
+      Field<DataType, CoordinateType> *me = this;
+
+      for(int xoffset=0; xoffset < hiresParallelisationSeparation; xoffset++) {
+#pragma omp parallel for default(none) shared(hires, me, pixel_volume_ratio, xoffset, hiresGridSize, hiresParallelisationSeparation)
+        for(int x=xoffset; x<hiresGridSize; x+=hiresParallelisationSeparation) {
+          for(int y=0; y<hiresGridSize; y++) {
+            for(int z=0; z<hiresGridSize; z++) {
+              Coordinate<int> coord {x,y,z};
+              auto index = hires.getGrid().getIndexFromCoordinateNoWrap(coord);
+              auto location = hires.getGrid().getCentroidFromCoordinate(coord);
+              me->deInterpolate(location, hires[index]/pixel_volume_ratio);
+            }
+          }
+        }
+      }
+    }
+
+
+
     //! Evaluates the field at the specified co-ordinate using interpolation.
     DataType evaluateInterpolated(Coordinate<CoordinateType> location) const {
       int x_p_0, y_p_0, z_p_0;

genetIC/src/simulation/field/multilevelfield.hpp

@@ -59,15 +59,14 @@ namespace fields {
     }
 
     //! Copy constructor
-    MultiLevelField(const MultiLevelField<DataType> &copy) :
-      std::enable_shared_from_this<MultiLevelField<DataType>>(), multiLevelContext(&(copy.getContext())) {
-
+    MultiLevelField(const MultiLevelField<DataType> &other) :
+      std::enable_shared_from_this<MultiLevelField<DataType>>(), multiLevelContext(&(other.getContext())) {
       for (size_t level = 0; level < multiLevelContext->getNumLevels(); level++) {
-        fieldsOnLevels.push_back(std::make_shared<Field<DataType, T>>(copy.getFieldForLevel(level)));
+        fieldsOnLevels.push_back(other.getFieldForLevel(level).copy());
       }
-      transferType = copy.transferType;
+      transferType = other.transferType;
 
-      isCovector = copy.isCovector;
+      isCovector = other.isCovector;
     }
 
     virtual ~MultiLevelField() {}
@@ -77,6 +76,11 @@ namespace fields {
       return transferType;
     }
 
+    //! Manually set the transfer function
+    void setForceTransferType(particle::species transfer_type) {
+      transferType = transfer_type;
+    }
+
     //! Returns a reference to the multi-level context associated to this multi-level field.
     virtual multilevelgrid::MultiLevelGrid<DataType> &getContext() const {
       return const_cast<multilevelgrid::MultiLevelGrid<DataType> &>(*multiLevelContext);
@@ -168,6 +172,11 @@ namespace fields {
       addScaled(other, 1.0);
     }
 
+    //! Subtracts the specified multi-level field from this one.
+    void operator-=(const MultiLevelField<DataType> &other) {
+      addScaled(other, -1.0);
+    }
+
     //! Divides the field by the specified ratio.
     void operator/=(DataType ratio) {
       using namespace tools::numerics;
@@ -178,6 +187,16 @@ namespace fields {
       }
     }
 
+    //! Multiplied the field by the specified ratio.
+    void operator*=(DataType ratio) {
+      using namespace tools::numerics;
+
+      for (size_t level = 0; level < getNumLevels(); level++) {
+        auto &data = getFieldForLevel(level).getDataVector();
+        data *= ratio;
+      }
+    }
+
     //! Flips the sign of the field.
     void reverse() {
       for (size_t level = 0; level < getNumLevels(); ++level) {
@@ -193,22 +212,32 @@ namespace fields {
     //! Add a scaled multilevel field to the current one
     void addScaled(const MultiLevelField &other, DataType scale) {
       assertContextConsistent();
-      assert(other.isFourierOnAllLevels());
       assert (isCompatible(other));
-      toFourier();
-
-      for (size_t level = 0; level < getNumLevels(); level++) {
-        if (hasFieldForLevel(level) && other.hasFieldForLevel(level)) {
+      if (other.isFourierOnAllLevels()) {
+        toFourier();
+
+        for (size_t level = 0; level < getNumLevels(); level++) {
+          if (hasFieldForLevel(level) && other.hasFieldForLevel(level)) {
+            Field<DataType> &fieldThis = getFieldForLevel(level);
+            const Field<DataType> &fieldOther = other.getFieldForLevel(level);
+            T kMin = fieldThis.getGrid().getFourierKmin();
+            fieldThis.forEachFourierCellInt([&fieldOther, kMin, scale]
+                                              (ComplexType currentVal, int kx, int ky, int kz) {
+              return currentVal + scale * fieldOther.getFourierCoefficient(kx, ky, kz);
+            });
+          }
+        }
+      } else if (other.isRealOnAllLevels()) {
+        toReal();
+        for (size_t level = 0; level < getNumLevels(); level++) {
           Field<DataType> &fieldThis = getFieldForLevel(level);
           const Field<DataType> &fieldOther = other.getFieldForLevel(level);
-          T kMin = fieldThis.getGrid().getFourierKmin();
-          fieldThis.forEachFourierCellInt([&fieldOther, kMin, scale]
-                                            (ComplexType currentVal, int kx, int ky, int kz) {
-            return currentVal + scale * fieldOther.getFourierCoefficient(kx, ky, kz);
-          });
+
+          fieldThis.addScaled(fieldOther, scale);
         }
+      } else {
+        throw std::runtime_error("Expected the other field to be in real/fourier space, not a mix of these.");
       }
-
     }
 
     //! \brief Copy across data from another field.
@@ -568,9 +597,26 @@ namespace fields {
       return *(this->fieldsOnLevels[i]);
     }
 
+    //! Combine all grids together
+    void combineGrids() {
+      auto filters = this->getFilters();
+      size_t nlevels = this->getContext().getNumLevels();
 
-  };
+      logging::entry() << "Combining information from different levels..." << std::endl;
+
+      for (size_t level = 1; level < nlevels; ++level) {
 
+        // remove the low-frequency information from this level
+        this->getFieldForLevel(level).applyFilter(
+          filters.getHighPassFilterForLevel(level));
+
+        // replace with the low-frequency information from the level below
+        this->getFieldForLevel(level).addFieldFromDifferentGridWithFilter(
+          getFieldForLevel(level - 1),
+          filters.getLowPassFilterForLevel(level - 1));
+      }
+   }
+  };
 
   /*! \class ConstraintField
     \brief Fields used for defining constraints.