DirectComputation.cpp 7.98 KB
Newer Older
1
// See LICENCE file at project root
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

#include <iostream>
#include <iomanip>

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>

#include  "ScalFmmConfig.h"
#include "Utils/FTic.hpp"
#include "Utils/FParameters.hpp"

#include "Files/FFmaGenericLoader.hpp"
#include "Kernels/P2P/FP2P.hpp"
17
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
18

19
#include "Utils/FParameterNames.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
20 21 22 23 24 25 26 27 28 29 30 31 32 33
//
/// \file  DirectComputation.cpp
//!
//! \brief DirectComputation: Driver to compute direct interaction between N particles for 1/r kernel.
//!
//! DirectComputation: Driver to compute direct interaction between N particles for 1/r kernel.
//! the particles are read from file given by -fin argument and potential, forces are stored in FMA format.
//!  <b> General arguments:</b>
//!     \param   -help (-h)      to see the parameters available in this driver
//!     \param   -fin name:  file name  to convert (with extension .fma (ascii) or bfma (binary).
//!                             Only our FMA (.bma, .bfma) is allowed "
//!     \param    -fout filenameOUT   output file  with extension (default output.bfma)
//!      \param   -verbose : print index x y z Q V fx fy fz
//!
34 35 36

// Simply create particles and try the kernels
int main(int argc, char ** argv){
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
  FHelpDescribeAndExit(argc, argv,
		       ">> This executable has to be used to compute  interaction either for periodic or non periodic system.\n"
		       ">> Example -fin filenameIN.{fma or bfma)     -fout filenameOUT{fma or bfma) \n"
		       ">> Default input file : Data/unitCubeXYZQ20k.fma\n"
		       ">> Only our FMA (.bma, .bfma) is allowed as input.\n"
		       ">> Output file  with extension (default output.bfma).",
		       FParameterDefinitions::InputFile, FParameterDefinitions::OutputFile,
		       FParameterDefinitions::EnabledVerbose);

  //////////////////////////////////////////////////////////////
  typedef double FReal;
  const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/unitCubeXYZQ20k.fma");
  const std::string filenameIn(FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options,  defaultFile.c_str()));
  const std::string filenameOut(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.bfma"));
  //
  FTic counter;

  // -----------------------------------------------------
  //  LOADER
  //  -----------------------------------------------------
  // ---------------------------------------------------------------------------------
  // Read  particles in the Octree
  // ---------------------------------------------------------------------------------
  std::cout << "Opening : " << filenameIn << "\n";
  //
  FFmaGenericLoader<FReal> loader(filenameIn);
  //
  FSize nbParticles = static_cast<int>(loader.getNumberOfParticles());
  std::cout << "Read " << nbParticles << " particles ..." << std::endl;
  double BoxWith=loader.getBoxWidth();
  FPoint<FReal> Centre(loader.getCenterOfBox().getX(), loader.getCenterOfBox().getY() , loader.getCenterOfBox().getZ());
  std::cout << "\tWidth : " <<BoxWith << " \t center x : " << loader.getCenterOfBox().getX()
	    << " y : " << loader.getCenterOfBox().getY() << " z : " << loader.getCenterOfBox().getZ() << std::endl;

  counter.tic();
  //
73 74
   FmaRWParticle<FReal, 4,8> *  particles = new FmaRWParticle<FReal, 4,8>[nbParticles]{};
 // memset(particles, 0, sizeof(FmaRWParticle<FReal, 4,8>) * nbParticles) ;
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
  //
  double totalCharge = 0.0;
  //
  //	int nbDataToRead = particles[0].getReadDataNumber();
  for(int idx = 0 ; idx<nbParticles ; ++idx){
    //
    loader.fillParticle(particles[idx].getPtrFirstData(), particles[idx].getReadDataNumber());
    //	loader.fillParticle(particles[idx].getPtrFirstData(), nbDataToRead);    // OK
    //  loader.fillParticle(particles[idx]); // OK
    //    std::cout << idx <<"  "<<  particles[idx].getPosition() << " "<<particles[idx].getPhysicalValue() << " "<<particles[idx].getPotential()
    //			<<"  " << particles[idx].getForces()[0]<<"  " <<particles[idx].getForces()[1]<<"  " <<particles[idx].getForces()[2]<<"  " <<std::endl;
    //
    totalCharge += particles[idx].getPhysicalValue() ;
  }

  counter.tac();

  std::cout << std::endl;
  std::cout << "Total Charge         = "<< totalCharge <<std::endl;
  std::cout << std::endl;

  std::cout << "Done  " << "(@ reading Particles  " << counter.elapsed() << " s)." << std::endl;
  //
  // ----------------------------------------------------------------------------------------------------------
  //                                   COMPUTATION
  // ----------------------------------------------------------------------------------------------------------
  // interaction kernel evaluator
  typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
  const MatrixKernelClass MatrixKernel;
  FReal denergy = 0.0;
  //
  //  computation
  //
  {
    printf("Compute :\n");
110
    counter.tic();
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
#pragma omp parallel shared(nbParticles, particles,denergy)
    {
#pragma omp for
      for(int idxTarget = 0 ; idxTarget < nbParticles ; ++idxTarget){
	//
	// compute with all other except itself
	//
	// Compute force and potential between  particles[idxTarget] and particles inside the box
	//
	for(int idxOther = 0; idxOther < nbParticles ; ++idxOther){
	  if( idxOther != idxTarget ){
	    FP2P::NonMutualParticles(particles[idxTarget].getPosition().getX(), particles[idxTarget].getPosition().getY(),
				     particles[idxTarget].getPosition().getZ(),particles[idxTarget].getPhysicalValue(),
				     &particles[idxTarget].setForces()[0],&particles[idxTarget].setForces()[1],
				     &particles[idxTarget].setForces()[2],particles[idxTarget].setPotential(),
				     particles[idxOther].getPosition().getX(), particles[idxOther].getPosition().getY(),
				     particles[idxOther].getPosition().getZ(),particles[idxOther].getPhysicalValue(),
				     &MatrixKernel);
	  }
	}
      } // end for
      // Compute the energy
#pragma omp  for reduction(+:denergy)
      for(int idx = 0 ; idx < nbParticles ; ++idx){
	denergy +=  particles[idx].getPotential()*(particles[idx].getPhysicalValue()) ;
      }
    } // end pragma parallel
138
    //
139
    denergy *= 0.5 ;
140 141
    counter.tac();
    //
142
    printf("Energy =   %.14e\n",denergy);
143
    //
144 145 146 147 148
    std::cout << "Done  " << "(@ Direct computation done = " << counter.elapsed() << " s)." << std::endl;
    std::cout << "\n"<< "END  "
	      << "-------------------------------------------------------------------------"
	      << std::endl << std::endl ;
  } // END
149 150 151 152 153 154 155

    //
    // ----------------------------------------------------------------
    //  Save  computation in binary format
    //
    //

156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
  std::cout << "Generate " << filenameOut <<"  for output file" << std::endl;
  //
  std::cout << " nbParticles: " << nbParticles <<"  " << sizeof(nbParticles) <<std::endl;
  std::cout << " denergy: " << denergy <<"  " << sizeof(denergy) <<std::endl;
  std::cout << " Box size: " << loader.getBoxWidth() << "  " << sizeof(loader.getBoxWidth())<<std::endl;
  //
  FFmaGenericWriter<FReal> writer(filenameOut) ;
  writer.writeHeader(Centre,BoxWith, nbParticles,*particles) ;
  writer.writeArrayOfParticles(particles, nbParticles);
  //
  // end generate
  // -----------------------------------------------------
  //
  if(FParameters::existParameter(argc, argv, FParameterDefinitions::EnabledVerbose.options)){
    denergy = 0 ;
    for(int idx = 0 ; idx < nbParticles ; ++idx){
      std::cout << ">> index " << idx << std::endl;
      std::cout << " x   " << particles[idx].getPosition().getX() << " y  " << particles[idx].getPosition().getY() << " z  " << particles[idx].getPosition().getZ() << std::endl;
      std::cout << " Q   " << particles[idx].getPhysicalValue()   << " V  " << particles[idx].getPotential() << std::endl;
      std::cout << " fx  " << particles[idx].getForces()[0]       << " fy " << particles[idx].getForces()[1]       << " fz " << particles[idx].getForces()[2] << std::endl;
      std::cout << "\n";
      denergy +=  particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
178
    }
179 180 181 182 183
  }
  std::cout << " ENERGY " << denergy << std::endl;
  //
  delete[] particles;
  return 0;
184 185 186 187
}