Commit 2e4a1b40 authored by adelmann's avatar adelmann 🎗
Browse files

remove unused stuff

parent 5d92f99e
......@@ -724,9 +724,6 @@ projects/mc4-new/src/Parser/Token.cpp -text
projects/mc4-new/src/Parser/Token.h -text
projects/mc4-new/src/Parser/TokenStream.cpp -text
projects/mc4-new/src/Parser/TokenStream.h -text
projects/mc4-new/src/PwrSpec/CMakeLists.txt -text
projects/mc4-new/src/PwrSpec/PwrSpec.cc -text
projects/mc4-new/src/PwrSpec/PwrSpec.hh -text
projects/mc4-new/src/SimData.hh -text
projects/mc4-new/src/Utilities/ArithmeticError.cpp -text
projects/mc4-new/src/Utilities/ArithmeticError.h -text
......
......@@ -56,7 +56,6 @@ endmacro()
add_subdirectory (Parser)
add_subdirectory (Utilities)
add_subdirectory (Initializer)
add_subdirectory (PwrSpec)
add_subdirectory (ChargedParticles)
add_subdirectory (FieldSolvers)
add_subdirectory (DataSink)
......
......@@ -9,7 +9,7 @@ template <class T, unsigned int Dim> class FieldSolver;
#include "FieldSolvers/FieldSolver.hh"
#include "Parser/FileStream.h"
#include "SimData.hh"
#include "PwrSpec/PwrSpec.hh"
#include "DataSink/DataSink.h"
#include <iostream>
......
set (_SRCS
)
include_directories (
${CMAKE_CURRENT_SOURCE_DIR}
)
add_sources(${_SRCS})
#include "PwrSpec.hh"
template <class T, unsigned int Dim>
PwrSpec<T,Dim>::PwrSpec(ChargedParticles<T,Dim> *beam, SimData<T,Dim> simData):
layout_m(&(beam->getFieldLayout())),
simData_m(simData)
{
doInit();
}
template <class T, unsigned int Dim>
void PwrSpec<T,Dim>::doInit() {
gDomain_m = layout_m->getDomain(); // global domain
lDomain_m = layout_m->getLocalNDIndex(); // local domain
kmax_m = nint(sqrt(3*(simData_m.ng_comp*simData_m.ng_comp/4))) + 1 ;
ng_m = simData_m.ng_comp;
spectra1D_m = (T *) malloc (kmax_m*sizeof(T));
Nk_m = (int *) malloc (kmax_m*sizeof(int));
for (int k=0; k<kmax_m; ++k) {
spectra1D_m[k] = 0.0;
Nk_m[k]=0;
}
}
template <class T, unsigned int Dim>
void PwrSpec<T,Dim>::calc1D(ChargedParticles<T,Dim> *univ)
{
unsigned int kk;
for (int i=0;i<kmax_m;i++) {
Nk_m[i]=0;
spectra1D_m[i] = 0.0;
}
INFOMSG("Sum psp=real( ... " << sum(univ->rho_m) << " kmax= " << kmax_m << endl;);
/*
Here scatter_+^%^$#% takes place
*/
NDIndex<Dim> loop;
// Octant 1
for (int i=gDomain_m[0].min(); i<(gDomain_m[0].max()+1)/2; ++i) {
loop[0]=Index(i,i);
for (int j=gDomain_m[1].min(); j<(gDomain_m[1].max()+1)/2; ++j) {
loop[1]=Index(j,j);
for (int k=gDomain_m[2].min(); k<(gDomain_m[2].max()+1)/2; ++k) {
loop[2]=Index(k,k);
if (lDomain_m.contains(loop)) {
kk=(int)nint(sqrt(i*i+j*j+k*k));
if (kk>kmax_m) {
INFOMSG("kk>kmax_m kk= " << kk << endl);
kk = kmax_m;
}
spectra1D_m[kk] += univ->rho_m.localElement(loop);
Nk_m[kk]++;
}
}
}
}
// Octant 2
for (int i=(gDomain_m[0].max()+1)/2; i<=gDomain_m[0].max(); ++i) {
loop[0]=Index(i,i);
for (int j=gDomain_m[1].min(); j<(gDomain_m[1].max()+1)/2; ++j) {
loop[1]=Index(j,j);
for (int k=gDomain_m[2].min(); k<(gDomain_m[2].max()+1)/2; ++k) {
loop[2]=Index(k,k);
if (lDomain_m.contains(loop)) {
kk=(int)nint(sqrt((i-ng_m)*(i-ng_m)+(j*j)+(k*k)));
if (kk>kmax_m) {
INFOMSG("kk>kmax_m kk= " << kk << endl);
kk = kmax_m;
}
spectra1D_m[kk] += univ->rho_m.localElement(loop);
Nk_m[kk]++;
}
}
}
}
// Octant 3
for (int i=gDomain_m[0].min(); i<(gDomain_m[0].max()+1)/2; ++i) {
loop[0]=Index(i,i);
for (int j=(gDomain_m[1].max()+1)/2; j<=gDomain_m[1].max(); ++j) {
loop[1]=Index(j,j);
for (int k=gDomain_m[2].min(); k<(gDomain_m[2].max()+1)/2; ++k) {
loop[2]=Index(k,k);
if (lDomain_m.contains(loop)) {
kk=(int)nint(sqrt((i*i)+(j-ng_m)*(j-ng_m)+(k*k)));
if (kk>kmax_m) {
INFOMSG("kk>kmax_m kk= " << kk << endl);
kk = kmax_m;
}
spectra1D_m[kk] += univ->rho_m.localElement(loop);
Nk_m[kk]++;
}
}
}
}
//Octant 4
for (int i=gDomain_m[0].min(); i<(gDomain_m[0].max()+1)/2; ++i) {
loop[0]=Index(i,i);
for (int j=gDomain_m[1].min(); j<(gDomain_m[1].max()+1)/2; ++j) {
loop[1]=Index(j,j);
for (int k=(gDomain_m[2].max()+1)/2; k<=gDomain_m[2].max(); ++k) {
loop[2]=Index(k,k);
if (lDomain_m.contains(loop)) {
kk=(int)nint(sqrt((i*i)+(j*j)+(k-ng_m)*(k-ng_m)));
if (kk>kmax_m) {
INFOMSG("kk>kmax_m kk= " << kk << endl);
kk = kmax_m;
}
spectra1D_m[kk] += univ->rho_m.localElement(loop);
Nk_m[kk]++;
}
}
}
}
//Octant 5
for (int i=gDomain_m[0].min(); i<(gDomain_m[0].max()+1)/2; ++i) {
loop[0]=Index(i,i);
for (int j=(gDomain_m[1].max()+1)/2; j<=gDomain_m[1].max(); ++j) {
loop[1]=Index(j,j);
for (int k=(gDomain_m[2].max()+1)/2; k<=gDomain_m[2].max(); ++k) {
loop[2]=Index(k,k);
if (lDomain_m.contains(loop)) {
kk=(int)nint(sqrt((i*i)+(j-ng_m)*(j-ng_m)+(k-ng_m)*(k-ng_m)));
if (kk>kmax_m) {
INFOMSG("kk>kmax_m kk= " << kk << endl);
kk = kmax_m;
}
spectra1D_m[kk]+=univ->rho_m.localElement(loop);
Nk_m[kk]++;
}
}
}
}
//Octant 6
for (int i=(gDomain_m[0].max()+1)/2; i<=gDomain_m[0].max(); ++i) {
loop[0]=Index(i,i);
for (int j=(gDomain_m[1].max()+1)/2; j<=gDomain_m[1].max(); ++j) {
loop[1]=Index(j,j);
for (int k=gDomain_m[2].min(); k<(gDomain_m[2].max()+1)/2; ++k) {
loop[2]=Index(k,k);
if (lDomain_m.contains(loop)) {
kk=(int)nint(sqrt((i-ng_m)*(i-ng_m)+(j-ng_m)*(j-ng_m)+(k*k)));
if (kk>kmax_m) {
INFOMSG("kk>kmax_m kk= " << kk << endl);
kk = kmax_m;
}
spectra1D_m[kk]+=univ->rho_m.localElement(loop);
Nk_m[kk]++;
}
}
}
}
//Octant 7
for (int i=(gDomain_m[0].max()+1)/2; i<=gDomain_m[0].max(); ++i) {
loop[0]=Index(i,i);
for (int j=gDomain_m[1].min(); j<(gDomain_m[1].max()+1)/2; ++j) {
loop[1]=Index(j,j);
for (int k=(gDomain_m[2].max()+1)/2; k<=gDomain_m[2].max(); ++k) {
loop[2]=Index(k,k);
if (lDomain_m.contains(loop)) {
kk=(int)nint(sqrt((i-ng_m)*(i-ng_m)+(j*j)+(k-ng_m)*(k-ng_m)));
if (kk>kmax_m) {
INFOMSG("kk>kmax_m kk= " << kk << endl);
kk = kmax_m;
}
spectra1D_m[kk]+=univ->rho_m.localElement(loop);
Nk_m[kk]++;
}
}
}
}
//Octant 8
for (int i=(gDomain_m[0].max()+1)/2; i<=gDomain_m[0].max(); ++i) {
loop[0]=Index(i,i);
for (int j=(gDomain_m[1].max()+1)/2; j<=gDomain_m[1].max(); ++j) {
loop[1]=Index(j,j);
for (int k=(gDomain_m[2].max()+1)/2; k<=gDomain_m[2].max(); ++k) {
loop[2]=Index(k,k);
if (lDomain_m.contains(loop)) {
kk=(int)nint(sqrt((i-ng_m)*(i-ng_m)+(j-ng_m)*(j-ng_m)+(k-ng_m)*(k-ng_m)));
if (kk>kmax_m) {
INFOMSG("kk>kmax_m kk= " << kk << endl);
kk = kmax_m;
}
spectra1D_m[kk]+=univ->rho_m.localElement(loop);
Nk_m[kk]++;
}
}
}
}
INFOMSG("Loops done" << endl);
reduce( &(Nk_m[0]), &(Nk_m[0]) + kmax_m , &(Nk_m[0]) ,OpAddAssign());
reduce( &(spectra1D_m[0]), &(spectra1D_m[0]) + kmax_m, &(spectra1D_m[0]) ,OpAddAssign());
}
template <class T, unsigned int Dim>
void PwrSpec<T,Dim>::save1DSpectra(string fn)
{
Inform* fdip = new Inform(NULL,fn.c_str(),Inform::OVERWRITE,0);
Inform& fdi = *fdip;
setInform(fdi);
setFormat(9,1,0);
T tpiL = 8.0*atan(1.0)/simData_m.rL; // 2 pi / rL
// Renormalize power spectrum to match mc2.
int sumNk = 0;
for (int i=0; i<kmax_m;i++) {
sumNk+=Nk_m[i];
spectra1D_m[i] /= 1.0*Nk_m[i];
spectra1D_m[i] = spectra1D_m[i]*std::pow((T)(1.0*simData_m.np),(T)3.0);
}
fdi << "# " ;
for (int i=0; i < kmax_m; i++)
fdi << "Nk[" << i << "]= " << Nk_m[i] << " ";
fdi << " sum Nk= " << sumNk << endl;
for (int i=0; i<simData_m.np/2;i++) {
spectra1D_m[i] = spectra1D_m[i+1];
fdi << (i+1)*tpiL << " " << spectra1D_m[i]*std::pow((T)(simData_m.rL/simData_m.ng_comp),(T)3.0) << " " << endl;
}
delete fdip;
}
#ifndef PWR_SPEC_H
#define PWR_SPEC_H
#include <ChargedParticles/ChargedParticles.hh>
#include <iostream>
#include <fstream>
#include <strstream>
// #include "AppTypes/dcomplex.h"
// #include "FFT/FFT.h"
template <class T, unsigned int Dim>
class PwrSpec {
typedef typename ChargedParticles<T,Dim>::Center_t Center_t;
typedef typename ChargedParticles<T,Dim>::Mesh_t Mesh_t;
typedef typename ChargedParticles<T,Dim>::FieldLayout_t FieldLayout_t;
typedef typename ChargedParticles<T,Dim>::Vector_t Vector_t;
typedef typename ChargedParticles<T,Dim>::IntrplCIC_t IntrplCIC_t;
T*spectra1D_m;
int *Nk_m;
FieldLayout_t *layout_m;
/// global domain for the various fields
NDIndex<Dim> gDomain_m;
/// local domain for the various fields
NDIndex<Dim> lDomain_m;
unsigned int ng_m;
unsigned int kmax_m;
SimData<T,Dim> simData_m;
/// fortrans nint function
inline T nint(T x)
{
return ceil(x + 0.5) - (fmod(x*0.5 + 0.25, 1.0) != 0);
}
void doInit();
public:
PwrSpec(ChargedParticles<T,Dim> *beam, SimData<T,Dim> simData);
~PwrSpec();
void calc1D(ChargedParticles<T,Dim> *b);
void save1DSpectra(string fn);
};
#include "PwrSpec.cc"
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment