CouplingCellPlotter.h

// Copyright (C) 2015 Technische Universitaet Muenchen
// This file is part of the Mamico project. For conditions of distribution
// and use, please see the copyright notice in Mamico's main folder, or at
// www5.in.tum.de/mamico
#ifndef _MOLECULARDYNAMICS_COUPLING_COUPLINGCELLPLOTTER_H_
#define _MOLECULARDYNAMICS_COUPLING_COUPLINGCELLPLOTTER_H_
 
#include "coupling/KineticEnergyController.h"
#include "coupling/cell-mappings/ComputeMassMapping.h"
#include "coupling/cell-mappings/ComputeMomentumMapping.h"
#include "coupling/cell-mappings/VTKMoleculePlotter.h"
#include "coupling/datastructures/LinkedCellContainer.h"
#include "coupling/interface/MDSolverInterface.h"
#include "tarch/la/Vector.h"
#include <fstream>
#include <sstream>
 
namespace coupling {
template <class LinkedCell, unsigned int dim> class CouplingCellPlotter;
}

template <class LinkedCell, unsigned int dim> class coupling::CouplingCellPlotter {
public:
  CouplingCellPlotter(unsigned int ID, std::string filename, unsigned int rank, unsigned int t, datastructures::LinkedCellContainer<LinkedCell, dim>& cells,
                      coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface) {
    plotMoleculeFile(ID, filename, rank, t, cells, mdSolverInterface);
    plotCouplingCellFile(ID, filename, rank, t, cells, mdSolverInterface);
  }

  ~CouplingCellPlotter() {}
 
private:
  void plotMoleculeFile(unsigned int ID, std::string filename, unsigned int rank, unsigned int t,
                        coupling::datastructures::LinkedCellContainer<LinkedCell, dim>& cells,
                        coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface) {
    // open files
    std::ofstream moleculeFile;
    // store positions and velocities of molecules
    std::stringstream moleculeVelocities;
    std::stringstream moleculePositions;
    std::stringstream moleculePotentials;
    // zeros appended in vtk file
    std::string appendFloatZeros("");
    switch (dim) {
    case 1:
      appendFloatZeros = " 0.0 0.0";
      break;
    case 2:
      appendFloatZeros = " 0.0";
      break;
    case 3:
      appendFloatZeros = "";
      break;
    default:
      break;
    }
    Writer4Molecules writer(moleculeVelocities, moleculePositions, moleculePotentials, appendFloatZeros, mdSolverInterface);
    cells.applyToLocalNonGhostCouplingCellsWithLinkedCells(writer);
    // open file and write info
    open(ID, filename, "Molecules", rank, t, moleculeFile);
    moleculeFile << "# vtk DataFile Version 2.0" << std::endl;
    moleculeFile << "Generated by MaMiCo (Philipp Neumann)" << std::endl;
    moleculeFile << "ASCII" << std::endl;
    moleculeFile << "DATASET UNSTRUCTURED_GRID" << std::endl;
    moleculeFile << "POINTS " << writer.getMoleculeCounter() << " float" << std::endl;
    moleculeFile << moleculePositions.str() << std::endl << std::endl;
    moleculeFile << "POINT_DATA " << writer.getMoleculeCounter() << std::endl;
    moleculeFile << "VECTORS velocity float" << std::endl;
    moleculeFile << moleculeVelocities.str() << std::endl;
    moleculeFile << "SCALARS potentialEnergy float 1" << std::endl;
    moleculeFile << "LOOKUP_TABLE DEFAULT" << std::endl;
    moleculeFile << moleculePotentials.str() << std::endl;
    moleculePositions.clear();
    moleculePositions.str("");
    moleculeVelocities.clear();
    moleculeVelocities.str("");
    moleculePotentials.clear();
    moleculePotentials.str("");
    moleculeFile.close();
  }

  void plotCouplingCellFile(unsigned int ID, std::string filename, unsigned int rank, unsigned int t,
                            coupling::datastructures::LinkedCellContainer<LinkedCell, dim>& cells,
                            coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface) {
    // stringstreams for all cell properties of interest
    std::stringstream microscopicMasses;
    microscopicMasses << "SCALARS microscopicMassBuffer float 1" << std::endl;
    microscopicMasses << "LOOKUP_TABLE DEFAULT" << std::endl;
    std::stringstream macroscopicMasses;
    macroscopicMasses << "SCALARS macroscopicMassBuffer float 1" << std::endl;
    macroscopicMasses << "LOOKUP_TABLE DEFAULT" << std::endl;
    std::stringstream microscopicMomenta;
    microscopicMomenta << "VECTORS microscopicMomentumBuffer float" << std::endl;
    std::stringstream macroscopicMomenta;
    macroscopicMomenta << "VECTORS macroscopicMomentumBuffer float" << std::endl;
    std::stringstream meanVelocities;
    meanVelocities << "VECTORS meanVelocity float" << std::endl;
    std::stringstream masses;
    masses << "SCALARS mass float 1" << std::endl;
    masses << "LOOKUP_TABLE DEFAULT" << std::endl;
    std::stringstream energies;
    energies << "SCALARS meanPotentialEnergyBuffer float 1" << std::endl;
    energies << "LOOKUP_TABLE DEFAULT" << std::endl;
    std::stringstream temperatures;
    temperatures << "SCALARS temperature float 1" << std::endl;
    temperatures << "LOOKUP_TABLE DEFAULT" << std::endl;
    // ofstream for file
    std::ofstream ofCell;
    // zeros appended in vtk file
    std::string appendFloatZeros("");
    switch (dim) {
    case 1:
      appendFloatZeros = " 0.0 0.0";
      break;
    case 2:
      appendFloatZeros = " 0.0";
      break;
    case 3:
      appendFloatZeros = "";
      break;
    default:
      break;
    }
 
    tarch::la::Vector<3, double> original = I11{{0, 0, 0}}.getCellMidPoint() - 0.5 * IDXS.getCouplingCellSize();
 
    // open file
    open(ID, filename, "CouplingCells", rank, t, ofCell);
    Writer4Cells writer(microscopicMasses, macroscopicMasses, microscopicMomenta, macroscopicMomenta, meanVelocities, masses, energies, temperatures,
                        appendFloatZeros, mdSolverInterface);
    cells.applyToAllLocalCouplingCellsWithLinkedCells(writer);
    ofCell << "# vtk DataFile Version 2.0" << std::endl;
    ofCell << "generated by MaMiCo (Philipp Neumann)" << std::endl;
    ofCell << "ASCII" << std::endl << std::endl;
    ofCell << "DATASET STRUCTURED_POINTS" << std::endl;
    ofCell << "DIMENSIONS " << I03::numberCellsInDomain << std::endl;
    ofCell << "ORIGIN ";
    for (unsigned int d = 0; d < 3; d++) {
      ofCell << original[d] << " ";
    }
    ofCell << std::endl;
    ofCell << "SPACING ";
    for (unsigned int d = 0; d < dim; d++) {
      ofCell << IDXS.getCouplingCellSize()[d] << " ";
    }
    for (unsigned int d = dim; d < 3; d++) {
      ofCell << "1.0 ";
    }
    ofCell << std::endl;
    ofCell << "POINT_DATA " << I02::linearNumberCellsInDomain << std::endl << std::endl;
    ofCell << microscopicMasses.str() << std::endl << std::endl;
    microscopicMasses.clear();
    microscopicMasses.str("");
    ofCell << macroscopicMasses.str() << std::endl << std::endl;
    macroscopicMasses.clear();
    macroscopicMasses.str("");
    ofCell << masses.str() << std::endl << std::endl;
    masses.clear();
    masses.str("");
    ofCell << energies.str() << std::endl << std::endl;
    energies.clear();
    energies.str("");
    ofCell << temperatures.str() << std::endl << std::endl;
    temperatures.clear();
    temperatures.str("");
    ofCell << microscopicMomenta.str() << std::endl << std::endl;
    microscopicMomenta.clear();
    microscopicMomenta.str("");
    ofCell << macroscopicMomenta.str() << std::endl << std::endl;
    macroscopicMomenta.clear();
    macroscopicMomenta.str("");
    ofCell << meanVelocities.str() << std::endl << std::endl;
    meanVelocities.clear();
    meanVelocities.str("");
    // close file
    ofCell.close();
  }

  void open(unsigned int ID, std::string filename, std::string fileType, unsigned int rank, unsigned int t, std::ofstream& of) {
    std::stringstream ss;
    ss << filename << "_" << fileType << "_" << ID << "_" << rank << "_" << t << ".vtk";
    of.open(ss.str().c_str());
    if (!of.is_open()) {
      std::cout << "ERROR coupling::CouplingCellPlotter: Could not open file " << ss.str() << "!" << std::endl;
      exit(EXIT_FAILURE);
    }
  }

  class Writer4Cells {
  public:
    Writer4Cells(std::stringstream& microscopicMasses, std::stringstream& macroscopicMasses, std::stringstream& microscopicMomenta,
                 std::stringstream& macroscopicMomenta, std::stringstream& meanVelocities, std::stringstream& masses, std::stringstream& energies,
                 std::stringstream& temperatures, const std::string& appendFloatZeros,
                 coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface)
        : _microscopicMasses(microscopicMasses), _macroscopicMasses(macroscopicMasses), _microscopicMomenta(microscopicMomenta),
          _macroscopicMomenta(macroscopicMomenta), _meanVelocities(meanVelocities), _masses(masses), _energies(energies), _temperatures(temperatures),
          _appendFloatZeros(appendFloatZeros), _computeMassMapping(mdSolverInterface), _computeMomentumMapping(mdSolverInterface),
          _kineticEnergyController(mdSolverInterface) {}

    ~Writer4Cells() {}

    void beginCellIteration() {}

    void endCellIteration() {}

    void apply(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, const I02& index) {
      // compute local quantities from molecules, if this is an inner cell
      double mass = 0.0;
      tarch::la::Vector<dim, double> meanVelocity(0.0);
      double temperature = 0.0;
 
      if (I10::contains(index)) {
        cell.iterateConstCells(_computeMassMapping);
        mass = _computeMassMapping.getMass();
        cell.iterateConstCells(_computeMomentumMapping);
        meanVelocity = _computeMomentumMapping.getMeanVelocity();
        temperature = _kineticEnergyController.computeTemperature(cell);
      }
 
      // write data to streams
      _microscopicMasses << cell.getMicroscopicMass() << std::endl;
      _macroscopicMasses << cell.getMacroscopicMass() << std::endl;
      _masses << mass << std::endl;
      _temperatures << temperature << std::endl;
      _energies << cell.getPotentialEnergy() << std::endl;
      for (unsigned int d = 0; d < dim; d++) {
        _microscopicMomenta << cell.getMicroscopicMomentum()[d] << " ";
        _macroscopicMomenta << cell.getMacroscopicMomentum()[d] << " ";
        _meanVelocities << meanVelocity[d] << " ";
      }
      _microscopicMomenta << _appendFloatZeros << std::endl;
      _macroscopicMomenta << _appendFloatZeros << std::endl;
      _meanVelocities << _appendFloatZeros << std::endl;
    }
 
  private:
    std::stringstream& _microscopicMasses;
    std::stringstream& _macroscopicMasses;
    std::stringstream& _microscopicMomenta;
    std::stringstream& _macroscopicMomenta;
    std::stringstream& _meanVelocities;
    std::stringstream& _masses;
    std::stringstream& _energies;
    std::stringstream& _temperatures;
    const std::string _appendFloatZeros;
    coupling::cellmappings::ComputeMassMapping<LinkedCell, dim> _computeMassMapping;
    coupling::cellmappings::ComputeMomentumMapping<LinkedCell, dim> _computeMomentumMapping;
    coupling::KineticEnergyController<LinkedCell, dim> _kineticEnergyController;
  };

  class Writer4Molecules {
  public:
    Writer4Molecules(std::stringstream& moleculeVelocities, std::stringstream& moleculePositions, std::stringstream& moleculePotentials,
                     const std::string& appendFloatZeros, coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface)
        : _moleculeVelocities(moleculeVelocities), _moleculePositions(moleculePositions), _moleculePotentials(moleculePotentials),
          _appendFloatZeros(appendFloatZeros), _mdSolverInterface(mdSolverInterface), _moleculeCounter(0) {}

    void beginCellIteration() {}

    void endCellIteration() {}

    void apply(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, const I02& index) {
      coupling::cellmappings::VTKMoleculePlotter<LinkedCell, dim> vtkMoleculePlotter(_moleculeVelocities, _moleculePositions, _moleculePotentials,
                                                                                     _appendFloatZeros, _mdSolverInterface);
      cell.iterateConstCells(vtkMoleculePlotter);
      _moleculeCounter += vtkMoleculePlotter.getParticleCounter();
    }

    unsigned int getMoleculeCounter() const { return _moleculeCounter; }
 
  private:
    std::stringstream& _moleculeVelocities;
    std::stringstream& _moleculePositions;
    std::stringstream& _moleculePotentials;
    const std::string _appendFloatZeros;
    coupling::interface::MDSolverInterface<LinkedCell, dim>* const _mdSolverInterface;
    unsigned int _moleculeCounter;
  };
};
#endif // _MOLECULARDYNAMICS_COUPLING_COUPLINGCELLPLOTTER_H_