Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Swapnil/cleanup development #102

Merged
merged 4 commits into from
Oct 2, 2024
Merged

Swapnil/cleanup development #102

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
387fc08
Cleanup
shankar-1729 Sep 30, 2024
2e3a8b3
more cleanup, enable compiling Emu with and without MPI
shankar-1729 Oct 1, 2024
45b9993
Just use AMREX_USE_MPI instead of creating another precompiler defini…
shankar-1729 Oct 1, 2024
4eb0463
cleanup
shankar-1729 Oct 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 6 additions & 15 deletions Source/Evolve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,34 +168,25 @@ void interpolate_rhs_from_mesh(FlavoredNeutrinoContainer& neutrinos_rhs, const M
yes_nulib, helperVarsReal_nulib, helperVarsInt_nulib);



//---------------------------

//The following commented loop can be used to print information about each particle in case debug is needed in future.
/*
const int lev = 0;

#ifdef _OPENMP
#pragma omp parallel
#endif
for (FNParIter pti(neutrinos_rhs, lev); pti.isValid(); ++pti)
{
const int np = pti.numParticles();
//printf("(Evolve.cpp) number of particles: np = %d \n", np);
FlavoredNeutrinoContainer::ParticleType* pstruct = &(pti.GetArrayOfStructs()[0]);

amrex::ParallelFor (np, [=] AMREX_GPU_DEVICE (int i) {
FlavoredNeutrinoContainer::ParticleType& p = pstruct[i];

//printf("(Inside Evolve.cpp)i =%d, Vphase = %g \n", i, p.rdata(PIdx::Vphase));
/*p.pos(0) = p.rdata(PIdx::x);
p.pos(1) = p.rdata(PIdx::y);
p.pos(2) = p.rdata(PIdx::z);

p.rdata(PIdx::x) = p.pos(0);
p.rdata(PIdx::y) = p.pos(1);
p.rdata(PIdx::z) = p.pos(2);*/

//printf("(Inside Evolve.cpp) Partile i = %d, Vphase = %g \n", i, p.rdata(PIdx::Vphase));
});
}
//--------------------------
*/



amrex::MeshToParticle(neutrinos_rhs, state, 0,
Expand Down
108 changes: 20 additions & 88 deletions Source/FlavoredNeutrinoContainerInit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -290,7 +290,6 @@ InitParticles(const TestParams* parms)

if(parms->IMFP_method == 1){
p.rdata(PIdx::Vphase) = dx[0]*dx[1]*dx[2]*4*MathConst::pi*(pow(p.rdata(PIdx::pupt)+parms->delta_E/2,3)-pow(p.rdata(PIdx::pupt)-parms->delta_E/2,3))/(3*ndirs_per_loc*parms->nppc[0]*parms->nppc[1]*parms->nppc[2]);
//printf("(Inside FlavoredNeutrinoContainerInit.cpp) Vphase = %g \n", p.rdata(PIdx::Vphase));
}

//=====================//
Expand Down Expand Up @@ -377,8 +376,6 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)
const auto dx = Geom(lev).CellSizeArray();
const auto plo = Geom(lev).ProbLoArray();
const auto& a_bounds = Geom(lev).ProbDomain();

//printf("plo = [%f, %f, %f] \n", plo[0], plo[1], plo[2]);

const int nlocs_per_cell = AMREX_D_TERM( parms->nppc[0],
*parms->nppc[1],
Expand All @@ -391,12 +388,8 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)

// determine the number of directions per location
int ndirs_per_loc = particle_data.size();
//amrex::Print() << "Using " << ndirs_per_loc << " directions." << std::endl;
const Real scale_fac = dx[0]*dx[1]*dx[2]/nlocs_per_cell;

//FIXME: We need to use outflow boundary condition, not periodic boundary condition. Put an assert(!periodic) here.


// Loop over multifabs //
#ifdef _OPENMP
#pragma omp parallel
Expand All @@ -415,16 +408,15 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)
const auto lo = amrex::lbound(tile_box);
const auto hi = amrex::ubound(tile_box);

//printf("tile_box = [%d, %d, %d] x [%d, %d, %d] \n", lo.x, lo.y, lo.z, hi.x, hi.y, hi.z);

Gpu::ManagedVector<unsigned int> counts(tile_box.numPts(), 0); //PODVector<int, ManagedArenaAllocator<int> > counts(n, 0)

unsigned int* pcount = counts.dataPtr();

Gpu::ManagedVector<unsigned int> offsets(tile_box.numPts());
unsigned int* poffset = offsets.dataPtr();

const int buffer = 0; //TODO: TODO: Set this appropriately.
//TODO: This can be used to emit particles not exactly at the boundary, but at n cells away from the boundary (n = buffer).
const int buffer = 0;

// Determine how many particles to add to the particle tile per cell
//This loop runs over all the particles in a given box.
Expand Down Expand Up @@ -479,8 +471,7 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)
}

if (!create_particle_this_cell) continue;
//printf("CREATE PARTRICLE AT: i = %d, j = %d, k = %d \n", i, j, k);


int ix = i - lo.x;
int iy = j - lo.y;
int iz = k - lo.z;
Expand All @@ -496,7 +487,6 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)
});

// Determine total number of particles to add to the particle tile
//Gpu::inclusive_scan(counts.begin(), counts.end(), offsets.begin()); //This sets the value of "offsets"
Gpu::exclusive_scan(counts.begin(), counts.end(), offsets.begin()); //This sets the value of "offsets"

int num_to_add = offsets[tile_box.numPts()-1] + counts[tile_box.numPts()-1];
Expand All @@ -521,9 +511,6 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)
offsets[i] += old_size;
}

//printf("num_to_add = %d, old_size = %lu, new_size = %lu \n", num_to_add, old_size, new_size);
//printf("Particle ID = %lu \n", ParticleType::NextID());

//Returns the next particle ID for this processor.
// Particle IDs start at 1 and are never reused. The pair, consisting of the ID and the CPU on which the particle is "born", is a globally unique identifier for a particle.
//The maximum of this value across all processors must be checkpointed and then restored on restart so that we don't reuse particle IDs.
Expand Down Expand Up @@ -559,11 +546,6 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)

get_position_unit_cell(r, parms->nppc, i_loc);

//Real x = plo[0] + (i + r[0])*dx[0];
//Real x = plo[0] for i_plus direction i.e VC, (y,z) remain same CC.
//Real x = plo[0] + (i+1)*dx[0] for i_minus direction i.e. VC, (y,z) remain same CC.
//Real y = plo[1] + (j + r[1])*dx[1];
//Real z = plo[2] + (k + r[2])*dx[2];
Real x, y, z;

bool create_particle_this_cell = false;
Expand Down Expand Up @@ -626,14 +608,10 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)
break;
}

//printf("x = %f, y = %f, z = %f \n", x, y, z);

if (!create_particle_this_cell) continue;
//printf("CREATE PARTRICLE AT: i = %d, j = %d, k = %d \n", i, j, k);

for(int i_direction=0; i_direction<ndirs_per_loc; i_direction++){
// Get the Particle data corresponding to our particle index in pidx
//const int pidx = poffset[cellid] - poffset[0] + i_loc*ndirs_per_loc + i_direction;
const int pidx = poffset[cellid] + i_loc*ndirs_per_loc + i_direction;
ParticleType& p = pstruct[pidx];

Expand Down Expand Up @@ -684,47 +662,41 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)
//p.rdata(PIdx::Vphase) = dx[0]*dx[1]*dx[2]*V_momentum;
const Real dt = current_dt;
const Real clight = PhysConst::c;
const Real pupx_ = p.rdata(PIdx::pupx); //TODO: Review this.
const Real pupy_ = p.rdata(PIdx::pupy); //TODO: Review this.
const Real pupz_ = p.rdata(PIdx::pupz); //TODO: Review this.
const Real pupt_ = p.rdata(PIdx::pupt); //TODO: Review this.
const Real pupx_ = p.rdata(PIdx::pupx);
const Real pupy_ = p.rdata(PIdx::pupy);
const Real pupz_ = p.rdata(PIdx::pupz);
const Real pupt_ = p.rdata(PIdx::pupt);

switch (DIRECTION)
{
//Create particles in +ve x direction at lower x boundary.
case BoundaryParticleCreationDirection::I_PLUS:
p.rdata(PIdx::Vphase) = dx[1]*dx[2]*clight*dt*V_momentum*std::abs(pupx_/pupt_);
//printf("(Inside FlavoredNeutrinoContainerInit.cpp:I_PLUS) Vphase = %g \n", p.rdata(PIdx::Vphase));
break;

//Create particles in -ve x direction at upper x boundary.
case BoundaryParticleCreationDirection::I_MINUS:
p.rdata(PIdx::Vphase) = dx[1]*dx[2]*clight*dt*V_momentum*std::abs(pupx_/pupt_);
//printf("(Inside FlavoredNeutrinoContainerInit.cpp:I_MINUS) Vphase = %g \n", p.rdata(PIdx::Vphase));
break;

//Create particles in +ve y direction at lower y boundary.
case BoundaryParticleCreationDirection::J_PLUS:
p.rdata(PIdx::Vphase) = dx[0]*dx[2]*clight*dt*V_momentum*std::abs(pupy_/pupt_);
//printf("(Inside FlavoredNeutrinoContainerInit.cpp:J_PLUS) Vphase = %g \n", p.rdata(PIdx::Vphase));
break;

//Create particles in -ve y direction at upper y boundary.
case BoundaryParticleCreationDirection::J_MINUS:
p.rdata(PIdx::Vphase) = dx[0]*dx[2]*clight*dt*V_momentum*std::abs(pupy_/pupt_);
//printf("(Inside FlavoredNeutrinoContainerInit.cpp:J_MINUS) Vphase = %g \n", p.rdata(PIdx::Vphase));
break;

//Create particles in +ve z direction at lower z boundary.
case BoundaryParticleCreationDirection::K_PLUS:
p.rdata(PIdx::Vphase) = dx[0]*dx[1]*clight*dt*V_momentum*std::abs(pupz_/pupt_);
//printf("(Inside FlavoredNeutrinoContainerInit.cpp:K_PLUS) Vphase = %g \n", p.rdata(PIdx::Vphase));
break;

//Create particles in -ve z direction at upper z boundary.
case BoundaryParticleCreationDirection::K_MINUS:
p.rdata(PIdx::Vphase) = dx[0]*dx[1]*clight*dt*V_momentum*std::abs(pupz_/pupt_);
//printf("(Inside FlavoredNeutrinoContainerInit.cpp:K_MINUS) Vphase = %g \n", p.rdata(PIdx::Vphase));
break;

default:
Expand All @@ -734,66 +706,26 @@ CreateParticlesAtBoundary(const TestParams* parms, const Real current_dt)
}
}

//TODO: Do not apply perturbation in this case.
//=====================//
// Apply Perturbations //
//=====================//
if(parms->perturbation_type == 0){
// random perturbations to the off-diagonals
/*Real rand;
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N01_Re) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N11_Re ));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N01_Im) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N11_Re ));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N01_Rebar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N11_Rebar));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N01_Imbar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N11_Rebar));*/
p.rdata(PIdx::N01_Re) = - p.rdata(PIdx::N11_Re);
p.rdata(PIdx::N01_Im) = - p.rdata(PIdx::N11_Re);
p.rdata(PIdx::N01_Rebar) = - p.rdata(PIdx::N11_Rebar);
p.rdata(PIdx::N01_Imbar) = - p.rdata(PIdx::N11_Rebar);
//Set off-diagonal terms to zero //TODO: This should be reviewed once.
p.rdata(PIdx::N01_Re) = 0.0;
p.rdata(PIdx::N01_Im) = 0.0;
p.rdata(PIdx::N01_Rebar) = 0.0;
p.rdata(PIdx::N01_Imbar) = 0.0;
#if NUM_FLAVORS==3
/*symmetric_uniform(&rand, engine);
p.rdata(PIdx::N02_Re) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N22_Re ));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N02_Im) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N22_Re ));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N12_Re) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N11_Re ) - p.rdata(PIdx::N22_Re ));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N12_Im) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N11_Re ) - p.rdata(PIdx::N22_Re ));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N02_Rebar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N22_Rebar));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N02_Imbar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N22_Rebar));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N12_Rebar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N11_Rebar) - p.rdata(PIdx::N22_Rebar));
symmetric_uniform(&rand, engine);
p.rdata(PIdx::N12_Imbar) = parms->perturbation_amplitude*rand * (p.rdata(PIdx::N11_Rebar) - p.rdata(PIdx::N22_Rebar));*/
p.rdata(PIdx::N02_Re) = - p.rdata(PIdx::N22_Re);
p.rdata(PIdx::N02_Im) = - p.rdata(PIdx::N22_Re);
p.rdata(PIdx::N12_Re) = - p.rdata(PIdx::N22_Re);
p.rdata(PIdx::N12_Im) = - p.rdata(PIdx::N22_Re);
p.rdata(PIdx::N02_Rebar) = - p.rdata(PIdx::N22_Rebar);
p.rdata(PIdx::N02_Imbar) = - p.rdata(PIdx::N22_Rebar);
p.rdata(PIdx::N12_Rebar) = - p.rdata(PIdx::N22_Rebar);
p.rdata(PIdx::N12_Imbar) = - p.rdata(PIdx::N22_Rebar);
p.rdata(PIdx::N02_Re) = 0.0;
p.rdata(PIdx::N02_Im) = 0.0;
p.rdata(PIdx::N12_Re) = 0.0;
p.rdata(PIdx::N12_Im) = 0.0;
p.rdata(PIdx::N02_Rebar) = 0.0;
p.rdata(PIdx::N02_Imbar) = 0.0;
p.rdata(PIdx::N12_Rebar) = 0.0;
p.rdata(PIdx::N12_Imbar) = 0.0;
#endif
}
if(parms->perturbation_type == 1){
// Perturb real part of e-mu component only sinusoidally in z
/*Real nu_k = (2.*M_PI) / parms->perturbation_wavelength_cm;
p.rdata(PIdx::N01_Re) = parms->perturbation_amplitude*sin(nu_k*p.pos(2)) * (p.rdata(PIdx::N00_Re ) - p.rdata(PIdx::N11_Re ));
p.rdata(PIdx::N01_Rebar) = parms->perturbation_amplitude*sin(nu_k*p.pos(2)) * (p.rdata(PIdx::N00_Rebar) - p.rdata(PIdx::N11_Rebar));*/
p.rdata(PIdx::N01_Re) = - p.rdata(PIdx::N11_Re);
p.rdata(PIdx::N01_Rebar) = - p.rdata(PIdx::N11_Rebar);
}

} // loop over direction
} // loop over location
}); // loop over grid cells

//printf("Finished creating particles at boundary. \n");
} // loop over multifabs

} // CreateParticlesAtBoundary()
Expand Down
8 changes: 0 additions & 8 deletions Source/NuLibTableHelpers.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,14 +72,6 @@ get_interp_spots_nulib(const double x, const double y, const double z, double &d
idx[6] = NTABLES_NULIB*(ix + nrho*((iy-1) + ntemp*((iz-1) + nye*(idx_species + nspecies*idx_group))));
idx[7] = NTABLES_NULIB*((ix-1) + nrho*((iy-1) + ntemp*((iz-1) + nye*(idx_species + nspecies*idx_group))));

//idx[0] = NTABLES_NULIB*(ix + nrho*(iy + ntemp*iz));
//idx[1] = NTABLES_NULIB*((ix-1) + nrho*(iy + ntemp*iz));
//idx[2] = NTABLES_NULIB*(ix + nrho*((iy-1) + ntemp*iz));
//idx[3] = NTABLES_NULIB*(ix + nrho*(iy + ntemp*(iz-1)));
//idx[4] = NTABLES_NULIB*((ix-1) + nrho*((iy-1) + ntemp*iz));
//idx[5] = NTABLES_NULIB*((ix-1) + nrho*(iy + ntemp*(iz-1)));
//idx[6] = NTABLES_NULIB*(ix + nrho*((iy-1) + ntemp*(iz-1)));
//idx[7] = NTABLES_NULIB*((ix-1) + nrho*((iy-1) + ntemp*(iz-1)));

// set up aux vars for interpolation
delx = logrho_nulib[ix] - x;
Expand Down
44 changes: 4 additions & 40 deletions Source/ReadEosTable.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,7 @@

#include "EosTable.H"

// mini NoMPI
#define HAVE_CAPABILITY_MPI //FIXME: This should be defined only when USE_MPI = TRUE
#ifdef HAVE_CAPABILITY_MPI
#ifdef AMREX_USE_MPI
#include <mpi.h>
#define BCAST(buffer, size) MPI_Bcast(buffer, size, MPI_BYTE, my_reader_process, MPI_COMM_WORLD)
#else
Expand Down Expand Up @@ -54,17 +52,10 @@ namespace nuc_eos_private {
void ReadEosTable(const std::string nuceos_table_name) {
using namespace nuc_eos_private;

//std::string nuceos_table_name = "/home/sshanka/000_UTK_projects/Emu/Exec/SFHo.h5";
amrex::Print() << "(ReadEosTable.cpp) Using table: " << nuceos_table_name << std::endl;

//TODO:
int my_reader_process = 0; //reader_process;
/*if (my_reader_process < 0 || my_reader_process >= CCTK_nProcs(cctkGH))
{
CCTK_VWarn(CCTK_WARN_COMPLAIN, __LINE__, __FILE__, CCTK_THORNSTRING,
"Requested IO process %d out of range. Reverting to process 0.", my_reader_process);
my_reader_process = 0;
}*/

const int read_table_on_single_process = 1;
//const int doIO = !read_table_on_single_process || CCTK_MyProc(cctkGH) == my_reader_process; //TODO:
Expand Down Expand Up @@ -187,16 +178,11 @@ void ReadEosTable(const std::string nuceos_table_name) {
// free memory of temporary array
myManagedArena.deallocate(alltables_temp, nrho_ * ntemp_ * nye_ * NTABLES);

// convert units, convert logs to natural log
// convert logs to natural log
// The latter is great, because exp() is way faster than pow()
// pressure
//energy_shift_ = energy_shift_ * EPSGF; //Old code.
//energy_shift_ = energy_shift_; //Let's not convert units yet.


for(int i=0;i<nrho_;i++) {
//logrho[i] = log(pow(10.0,logrho[i]) * RHOGF);
// by using log(a^b*c) = b*log(a)+log(c)
//logrho[i] = logrho[i] * log(10.) + log(RHOGF); //Old code.
logrho[i] = logrho[i] * log(10.); //Let's not convert units yet. Only convert log_10(rho) to ln(rho).
}

Expand All @@ -213,42 +199,20 @@ void ReadEosTable(const std::string nuceos_table_name) {
assert(0);
}

// convert units //FIXME: We do not convert units yet. Just convert log10 to natural log.
//convert log10 to natural log.
for(int i=0;i<nrho_*ntemp_*nye_;i++) {

{ // pressure
int idx = 0 + NTABLES*i;
//alltables[idx] = alltables[idx] * log(10.0) + log(PRESSGF); //old code
alltables[idx] = alltables[idx] * log(10.0); //Let's not convert units yet.
}

{ // eps
int idx = 1 + NTABLES*i;
//alltables[idx] = alltables[idx] * log(10.0) + log(EPSGF); //old code
alltables[idx] = alltables[idx] * log(10.0);
epstable[i] = exp(alltables[idx]); //Let's not convert units yet.
}

/*{ // cs2
int idx = 4 + NTABLES*i;
alltables[idx] *= LENGTHGF*LENGTHGF/TIMEGF/TIMEGF;
}

{ // dedT
int idx = 5 + NTABLES*i;
alltables[idx] *= EPSGF;
}

{ // dpdrhoe
int idx = 6 + NTABLES*i;
alltables[idx] *= PRESSGF/RHOGF;
}

{ // dpderho
int idx = 7 + NTABLES*i;
alltables[idx] *= PRESSGF/EPSGF;
}*/

}

//allocate memory for helperVars
Expand Down
Loading