TimeInterp.H

Go to the documentation of this file.

 
#ifndef _TIMEINTERP_H_
#define _TIMEINTERP_H_
 
#include <vector>
#include <string>
#include <regex>
#include <algorithm>
#include <filesystem>
#include <type_traits>
#include "CartesianGrid.H"
#include "GizmoGas.H"
#include "Gas.H"
#include "GasData.H"
 
namespace criptic {
  namespace gas {
 
    template <typename T> class TimeInterp : public Gas {
 
    static_assert(std::is_base_of<Gas, T>::value,
      "Template class of TimeInterp must be derived from Gas");
 
      public:
      
      TimeInterp(const ParmParser& pp,
        const Geometry& geom_,
        const IdxType nBB,
        const int nGhost,
        std::vector<Real> TBB_ = {},
        std::vector<Real> WBB_ = {}) {
        
        // Save verbosity
        verbosity = 1;
        pp.query("verbosity", verbosity);
 
        std::string dirName, tempName;
        if (pp.query("gas.slice_dir", dirName) && 
            pp.query("gas.slice_filename", tempName)) {
          std::regex reg(tempName);
        
          // Find all the file names that match the template; these are files
          // in the directory dirName with names that are regex matches to
          // slice_filename
          for (auto const& file : std::filesystem::directory_iterator(dirName)) {
            if (!file.is_regular_file()) continue;          // Skip non-files
            std::string fname = file.path().filename().string(); // Extract file name
            if (std::regex_search(fname, reg))                    // Check for match
            //snapFiles.push_back(dirName + file.path().string());
            snapFiles.push_back(file.path().string());
          }
        };
 
 
        // Read the time slices available; this can be specified either by
        // just giving the time between slices and the total number
        // available, or by giving a list of all times available
        // Read in gas.times to get tSlice. If it is not supplied, continue
        if (!pp.query("gas.times", tSlice)) {
          Real dt;
          int nSlice;
          // Handle case where nSlice and dt are supplied. Snapshot times will
          // be an integer number times dt
          if (pp.query("gas.dt", dt) && pp.query("gas.nSlice", nSlice)) {
            if (nSlice < 2) {
            MPIUtil::haltRun("TimeInterp: need >= 2 slice times for gas.nSlice",
                errBadInputFile);
            }
            tSlice.resize(nSlice);
            for (int i=0; i < nSlice; i++) tSlice[i] = i*dt;
          } else {
            tSlice.resize(snapFiles.size());
            for (IdxType i=0; i < snapFiles.size(); i++) {
              // Find times of snapshot files individually.
              // NOTE: T::loadTime must be supplied in the definition of T
              tSlice[i] = T::loadTime(snapFiles[i]);
            }
          }
        }
        // Sort snapshot files and time slices
        snapFiles = sortSnapFiles(snapFiles, tSlice); 
        std::sort(tSlice.begin(), tSlice.end());
 
        // Supply an offset time to start the simulation at
        Real slice_t0 = 0;
        bool t0_supplied = pp.query("gas.slice_t0", slice_t0);
 
        if (slice_t0<0) {
          MPIUtil::haltRun("TimeInterp: gas.slice_t0 is negative", errBadInputFile);
        }
        
        if (t0_supplied) {
          // Ensure t0 is before the last snapshot
          if (slice_t0 > tSlice[tSlice.size()-1]) {
            MPIUtil::haltRun("TimeInterp: gas.slice_t0 exceeds the time of the final snapshot file",
              errBadInputFile);
          }
          // Aim: find the index of a file whose next file is immediately before t0
          IdxType firstIdx = -1;
          for (IdxType i=0; i<tSlice.size(); i++) {
            if (tSlice[i] < slice_t0) {
              firstIdx += 1; // Slice time is below t0; shift index
            } else {
              tSlice[i] += -slice_t0; // Shift time on slices after t0
            }
          }
          if (firstIdx > 0) { // Only delete files if the second file's time is less than t0
            // Delete times until before the last file before t0
            tSlice.erase(begin(tSlice), begin(tSlice) + firstIdx);
            snapFiles.erase(begin(snapFiles), begin(snapFiles) + firstIdx);
          }
          // Shift the first time. This may be negative, but this is necessary
          // for interpolation to work at the first pair of snapshots
          tSlice[0] += -slice_t0;
        };
        
        // Make sure the first time slice is time zero, and that times are
        // monotonically increasing
        if (!t0_supplied && tSlice[0] != 0.0)
        MPIUtil::haltRun("TimeInterp: first time slice must be at t = 0",
        errBadInputFile);
        for (IdxType i = 0; i < tSlice.size()-1; i++) {
        if (tSlice[i+1] <= tSlice[i])
        MPIUtil::haltRun("TimeInterp: slice times must be "
          "monotonically increasing",
          errBadInputFile);
        }
 
        // Make sure initial time step is less than time of first slice; it
        // is an error if not
        Real dt = 0.0;
        pp.query("dt_init", dt);
        if (dt >= tSlice[1])
        MPIUtil::haltRun("TimeInterp: initial simulation "
        "time step must be < time of second slice",
        errBadInputFile);
        //Store maximum smoothing length and leaf size
        if (!pp.query("gas.tree.maxH", hLim)) {hLim = 0;}
        if (!pp.query("gas.tree.leafSize", leafSize)) {leafSize = 32;}
 
        // Store uniform radiation field data
        if (TBB_.size() != WBB_.size()) {
          MPIUtil::haltRun("TimeInterp: Radiation field temperatures and "
          "dilution factors must be vectors of the same size",
          errBadInputFile);
        }
        if (TBB_.size()>0) {
          // Uniform radiation field data is supplied
          uniformRadField = true;
        } else { 
          // Radiation fields are not supplied: either no fields are necessary
          // or they are defined in T::gasData
          uniformRadField = false;
        }
        TBB = TBB_;
        WBB = WBB_;
 
        // Initialise old and new gas data objects as empty; they will be filled
        // using T::loadData during TimeInterp<T>::updateState
        gasOld = nullptr;
        gasNew = nullptr;
        // Flag that no data have been loaded yet
        dataLoaded = false;
      }
 
      
      virtual ~TimeInterp() {
        delete gasOld;
        delete gasNew;
      }
      
      virtual GasData gasData(const RealVec &x,
                  const Real t) const override {
        Real w = (t - tSlice[tPtr-1]) / (tSlice[tPtr] - tSlice[tPtr-1]);
        GasData gdOld = gasOld->gasData(x, t);
        GasData gdNew = gasNew->gasData(x, t);
        GasData gd_interp = (1 - w) * gdOld + w * gdNew;
        if (uniformRadField) {
          gd_interp.TBB = TBB;
          gd_interp.WBB = WBB;
        }
        return gd_interp;
      }
 
      virtual Real dxGhost() const override {
        return gasOld->dxGhost();
      }
 
      void updateState(const Real t,
                   Real& tNext) {
        // Safety check
        if (t >= tSlice.back())
        MPIUtil::haltRun("TimeInterp: simulation time has "
        "exceeded time of last time slice",
        errBadInputFile);
 
        // Handle special case of the first time this method is called, at
        // which point no data are loaded and tPtr does not yet point to the
        // right location.
        if (!dataLoaded) {
 
        // Move pointer so tSlice[tPtr] <= t < tSlice[tPtr+1]
        tPtr = 0;
        while (tSlice[tPtr+1] < t) tPtr++;
 
        // Load old and new data. Boundary conditions on cartesian snapshots 
        // must be set inside of the class declaration itself and not here
        if (verbosity > 1 && MPIUtil::IOProc) {
          std::cout << "TimeInterp: loading file " << snapFiles[tPtr] <<
           " at t = " << tSlice[tPtr] << std::endl;
        }
 
        gasOld = T::loadData(snapFiles[tPtr], hLim, leafSize);
 
        if (verbosity > 1 && MPIUtil::IOProc) {
          std::cout << "TimeInterp: loading file " << snapFiles[tPtr+1] <<
           " at t = " << tSlice[tPtr+1] << std::endl;
        }
        tPtr++;
        gasNew = T::loadData(snapFiles[tPtr], hLim, leafSize);       
        // Flag that data are now loaded
        dataLoaded = true;
 
        // Return
        return;
        }
 
        // General case on second and subsequent calls to this method: check
        // if we have hit the time when we need to load a new snapshot, or
        // if we need to reduce the time step to hit the next load time
        if (t == tSlice[tPtr]) {
 
        std::swap(gasOld, gasNew);    // std::swap old and new data
        tPtr++; // Increment pointer
        if (verbosity > 1 && MPIUtil::IOProc) {
          std::cout << "TimeInterp: loading file " << snapFiles[tPtr] <<
           " at t = " << tSlice[tPtr] << std::endl;
        }
 
        // Replace gasNew with new data
        delete gasNew;
        gasNew = T::loadData(snapFiles[tPtr], hLim, leafSize);
        } else if (tNext > tSlice[tPtr]) {
        tNext = tSlice[tPtr];    // Throttle time step
        }  
      }
 
      static std::vector<std::string> sortSnapFiles(std::vector<std::string> files, 
        std::vector<Real> times) {
 
        std::vector<std::pair<double, IdxType>> indexedTimes;
        for (IdxType i=0; i<times.size(); i++) {
            indexedTimes.emplace_back(times[i], i);
        }
        std::sort(indexedTimes.begin(), indexedTimes.end());
        std::vector<std::string> sortedFiles(files.size());
        for (IdxType i=0; i<files.size(); i++) {
            sortedFiles[i] = files[indexedTimes[i].second];
        }  
        return sortedFiles;
      }
      
      protected:
      // Data
      std::vector<Real> tSlice; 
      std::vector<std::string> snapFiles; 
      IdxType tPtr; 
      Real hLim;    
      int leafSize;    
      bool uniformRadField; 
      std::vector<Real> TBB; 
      std::vector<Real> WBB; 
      private:
      // Data
      int verbosity;   
      T * gasOld;      
      T * gasNew;      
      bool dataLoaded; 
    };
  }
}
 
#endif
// _TIMEINTERP_H_