Splicing multilevel

genetIC/src/ic.hpp

@@ -1659,7 +1659,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(size_t newSeed, T accuracy) {
     initialiseRandomComponentIfUninitialised();
     if(outputFields.size()>1)
       throw std::runtime_error("Splicing is not yet implemented for the case of multiple transfer functions");
@@ -1680,14 +1680,16 @@ class ICGenerator {
     newGenerator.draw();
     logging::entry() << "Finished constructing new random field. Beginning splice operation." << 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));
-      splicedFieldThisLevel.toFourier();
-      originalFieldThisLevel = std::move(splicedFieldThisLevel);
-    }
+    modifications::splice(newField, *(outputFields[0]), accuracy);
+
+    // 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);
+    // }
   }
 
   //! Reverses the sign of the low-k modes.

genetIC/src/simulation/field/multilevelfield.hpp

@@ -168,6 +168,11 @@ namespace fields {
       addScaled(other, 1.0);
     }
 
+    //! Subtracts the specified multi-level field to 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 +183,15 @@ namespace fields {
       }
     }
 
+    void operator*=(DataType val) {
+      using namespace tools::numerics;
+
+      for (size_t level = 0; level < getNumLevels(); level++) {
+        auto &data = getFieldForLevel(level).getDataVector();
+        data *= val;
+      }
+    }
+
     //! Flips the sign of the field.
     void reverse() {
       for (size_t level = 0; level < getNumLevels(); ++level) {
@@ -193,19 +207,17 @@ namespace fields {
     //! Add a scaled multilevel field to the current one
     void addScaled(const MultiLevelField &other, DataType scale) {
       assertContextConsistent();
-      assert(other.isFourierOnAllLevels());
+      if (other.isFourierOnAllLevels()) toFourier();
+      else if (other.isRealOnAllLevels()) toReal();
+      else throw std::runtime_error("Incompatible field types");
+
       assert (isCompatible(other));
-      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);
-          });
+          fieldThis.addScaled(fieldOther, scale);
         }
       }
 
@@ -569,6 +581,14 @@ namespace fields {
     }
 
 
+    //! Returns a constant reference to the field on the specified grid, assuming it has been populated.
+    const Field<DataType, T> &getFieldForLevel(size_t i) const override {
+      this->assertContextConsistent();
+      assert(i < this->fieldsOnLevels.size());
+
+      return *(this->fieldsOnLevels[i]);
+    }
+
   };