Skip to content

Commit

Permalink
Merge branch 'flory_huggins' of github.com:AMReX-FHD/FHDeX into flory…
Browse files Browse the repository at this point in the history 
…_huggins
  • Loading branch information
@jbbel
jbbel committed Nov 13, 2023
2 parents 481f5bd + 46f74a3 commit 7bebd68
Show file tree
Hide file tree
Showing 683 changed files with 5,814,358 additions and 2,993 deletions.
215 changes: 99 additions & 116 deletions exec/DSMC/DsmcCollide.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,18 @@ int FhdParticleContainer::getSpeciesIndex(int species1, int species2)
}
}

AMREX_GPU_HOST_DEVICE AMREX_INLINE
void getSpeciesIndexRet(int species1, int species2, int* ret)
{
if(species1<species2)
{
ret[0] = species2+nspecies*species1;
} else
{
ret[0] = species1+nspecies*species2;
}
}

void FhdParticleContainer::InitCollisionCells()
{
BL_PROFILE_VAR("InitCollisionCells()",InitCollisionCells);
Expand Down Expand Up @@ -49,7 +61,8 @@ void FhdParticleContainer::InitCollisionCells()
int ij_spec;
for (int i_spec=0; i_spec<nspecies; i_spec++) {
for (int j_spec = i_spec; j_spec < nspecies; j_spec++) {
ij_spec = getSpeciesIndex(i_spec,j_spec);
//ij_spec = getSpeciesIndex(i_spec,j_spec);
getSpeciesIndexRet(i_spec,j_spec, &ij_spec);
arrselect(i,j,k,ij_spec) = 0.0;
}
}
Expand All @@ -73,8 +86,21 @@ void FhdParticleContainer::CalcSelections(Real dt)

const Array4<Real> & arrvrmax = mfvrmax.array(mfi);
const Array4<Real> & arrselect = mfselect.array(mfi);

amrex::ParallelFor(tile_box,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {

auto inds = m_bins.permutationPtr();
auto offs = m_bins.offsetsPtr();

Real ocollisionCellVolTmp = ocollisionCellVol;
Real particle_neff_tmp = particle_neff;

Real csx[MAX_SPECIES*MAX_SPECIES];

for(int i=0;i<(nspecies*nspecies);i++)
{
csx[i] = interproperties[i].csx;
}

amrex::ParallelForRNG(tile_box,[=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept {
int ij_spec;
long np_i, np_j;

Expand All @@ -88,21 +114,20 @@ void FhdParticleContainer::CalcSelections(Real dt)
for (int i_spec = 0; i_spec<nspecies; i_spec++)
{
for (int j_spec = i_spec; j_spec < nspecies; j_spec++) {
ij_spec = getSpeciesIndex(i_spec,j_spec);
np_i = m_cell_vectors[i_spec][grid_id][imap].size();
np_j = m_cell_vectors[j_spec][grid_id][imap].size();

//ij_spec = getSpeciesIndex(i_spec,j_spec);
getSpeciesIndexRet(i_spec,j_spec,&ij_spec);
//np_i = m_cell_vectors[i_spec][grid_id][imap].size();
//np_j = m_cell_vectors[j_spec][grid_id][imap].size();
np_i = getBinSize(offs,iv,i_spec,tile_box);
np_j = getBinSize(offs,iv,j_spec,tile_box);
vrmax = arrvrmax(i,j,k,ij_spec);
crossSection = interproperties[ij_spec].csx;
crossSection = csx[ij_spec];
//crossSection = 0;
if(i_spec==j_spec) {np_j = np_i-1;}
NSel = particle_neff*np_i*np_j*crossSection*vrmax*ocollisionCellVol*dt;
NSel = particle_neff_tmp*np_i*np_j*crossSection*vrmax*ocollisionCellVolTmp*dt;
if(i_spec==j_spec) {NSel = NSel*0.5;}
arrselect(i,j,k,ij_spec) = std::floor(NSel + amrex::Random());
arrselect(i,j,k,ij_spec) = std::floor(NSel + amrex::Random(engine));

if(i==1)
{
//Print() << "spec " << i_spec << " and " << j_spec << " selecting " << arrselect(i,j,k,ij_spec) << endl;
}
}
}
});
Expand All @@ -119,78 +144,81 @@ void FhdParticleContainer::CollideParticles(Real dt)
const int grid_id = mfi.index();
const int tile_id = mfi.LocalTileIndex();
auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)];
auto& particles = particle_tile.GetArrayOfStructs();
auto& aos = particle_tile.GetArrayOfStructs();
ParticleType* particles = aos().dataPtr();
//const long np = particle_tile.numParticles();
//auto& particles = particle_tile.GetArrayOfStructs();

const Array4<Real> & arrvrmax = mfvrmax.array(mfi);
const Array4<Real> & arrselect = mfselect.array(mfi);

auto inds = m_bins.permutationPtr();
auto offs = m_bins.offsetsPtr();


//const long np = particles.numParticles();
//amrex::ParallelForRNG(tile_box,
// [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept {
IntVect smallEnd = tile_box.smallEnd();
IntVect bigEnd = tile_box.bigEnd();

for (int i = smallEnd[0]; i <= bigEnd[0]; i++) {
for (int j = smallEnd[1]; j <= bigEnd[1]; j++) {
for (int k = smallEnd[2]; k <= bigEnd[2]; k++) {

dsmcSpecies* propertiesTmp = properties.data();

Real mass[MAX_SPECIES];

for(int i=0;i<(nspecies);i++)
{
mass[i] = properties[i].mass;
}

// for (int i = smallEnd[0]; i <= bigEnd[0]; i++) {
// for (int j = smallEnd[1]; j <= bigEnd[1]; j++) {
// for (int k = smallEnd[2]; k <= bigEnd[2]; k++) {
amrex::ParallelForRNG(tile_box,[=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept {
int totalCol = 0;
const IntVect& iv = {i,j,k};
long imap = tile_box.index(iv);

Real NSel[nspecies*nspecies], totalSel, selrun;
int np[nspecies];
int pindxi, pindxj;
Real NSel[MAX_SPECIES*MAX_SPECIES], totalSel, selrun;
int np[MAX_SPECIES];
int pindxi, pindxj, bini, binj;
int ij_spec;


for(int kk=0;kk<4;kk++)
{
// if(i==1)
//{
//NSel[kk] = 0;
//}
}

RealVect eij, vreij;
RealVect vi, vj, vij;
RealVect vijpost, boost;
Real massi, massj, massij;
Real vrmag, vrmax, vreijmag;

unsigned int* specLists[MAX_SPECIES];

for(int i=0;i<nspecies;i++)
{
specLists[i] = getCellList(inds,offs,iv,i,tile_box);
}

totalSel = 0;
// for (int i_spec = nspecies-1; i_spec>=0; i_spec--)
for(int i_spec = 0; i_spec<nspecies; i_spec++)
{
np[i_spec] = m_cell_vectors[i_spec][grid_id][imap].size();
//np[i_spec] = m_cell_vectors[i_spec][grid_id][imap].size();
np[i_spec] = getBinSize(offs,iv,i_spec,tile_box);
for (int j_spec = i_spec; j_spec < nspecies; j_spec++)
{
// int j_spec = i_spec;
ij_spec = getSpeciesIndex(i_spec,j_spec);
//ij_spec = getSpeciesIndex(i_spec,j_spec);
getSpeciesIndexRet(i_spec,j_spec, &ij_spec);
NSel[ij_spec] = (int)arrselect(i,j,k,ij_spec);
totalSel += NSel[ij_spec];
if(i==0)
{
//Print() << "cell " << i << " ij " << ij_spec << ": " << NSel[ij_spec] << endl;
}
if(i==127)
{
//Print() << "cell " << i << " ij " << ij_spec << ": " << NSel[ij_spec] << endl;
}

}

}
if(i==1)
{
//Print() << "total: " << totalSel << endl;
}


int totalSelect2 = totalSel;
int speci, specj, specij;


while (totalSel>0)
{
Real RR = amrex::Random();
Real RR = amrex::Random(engine);
bool spec_select = false;
selrun = 0;
speci = -1; specj = -1; specij = -1;
Expand All @@ -202,12 +230,9 @@ void FhdParticleContainer::CollideParticles(Real dt)
for(int j_spec=i_spec;j_spec<nspecies;j_spec++)
{
//int j_spec = i_spec;
ij_spec = getSpeciesIndex(i_spec,j_spec);
getSpeciesIndexRet(i_spec,j_spec,&ij_spec);
selrun += NSel[ij_spec];
if(i==1 && totalSelect2 == totalSel)
{
//Print() << "selrun: " << ij_spec << ": " << NSel[ij_spec] << endl;
}

if(selrun/totalSel>RR && !spec_select)
{
spec_select = true;
Expand All @@ -219,56 +244,21 @@ void FhdParticleContainer::CollideParticles(Real dt)
}
}
}


// Real specTest;
// int ci=0, cj=0;
// Real currentMax = 0;
// //for (int i_spec = nspecies-1; i_spec>=0; i_spec--)
// for(int i_spec = 0; i_spec<nspecies; i_spec++)
// {
// for(int j_spec=i_spec;j_spec<nspecies;j_spec++)
// {
// ij_spec = getSpeciesIndex(i_spec,j_spec);
// specTest = (NSel[ij_spec]/totalSel)*amrex::Random();
// //Print() << "SpecTest: " << specTest << endl;
// if(specTest > currentMax){ci = i_spec;cj = j_spec; currentMax = specTest;}
// }
// }
// //Print() << "ci, cj: " << ci << ", " << cj << endl;
// specij = getSpeciesIndex(ci,cj);
// speci = ci; specj = cj;
// NSel[ij_spec]--;


// Real specTest;
// int breaktest = 0;
// int ci=0, cj=0;
// Real currentMax = 0;
// for (int i_spec = nspecies-1; i_spec>=0; i_spec--)
// //for(int i_spec = 0; i_spec<nspecies; i_spec++)
// {
// for(int j_spec=i_spec;j_spec<nspecies;j_spec++)
// {
// ij_spec = getSpeciesIndex(i_spec,j_spec);
// if(NSel[ij_spec] > 0 && breaktest == 0){ci = i_spec;cj = j_spec; breaktest = 1;}
// }
// }
// //Print() << "ci, cj: " << ci << ", " << cj << endl;
// specij = getSpeciesIndex(ci,cj);
// speci = ci; specj = cj;
// NSel[ij_spec]--;


totalSel--;
massi = properties[speci].mass;
massj = properties[specj].mass;
massij = properties[speci].mass + properties[specj].mass;
massi = mass[speci];
massj = mass[specj];
massij = mass[speci] + mass[specj];
vrmax = arrvrmax(i,j,k,specij);
pindxi = (int)floor(amrex::Random()*m_cell_vectors[speci][grid_id][imap].size());
pindxj = (int)floor(amrex::Random()*m_cell_vectors[specj][grid_id][imap].size());
pindxi = m_cell_vectors[speci][grid_id][imap][pindxi];
pindxj = m_cell_vectors[specj][grid_id][imap][pindxj];
// pindxi = (int)floor(amrex::Random()*m_cell_vectors[speci][grid_id][imap].size());
// pindxj = (int)floor(amrex::Random()*m_cell_vectors[specj][grid_id][imap].size());
pindxi = (int)floor(amrex::Random(engine)*getBinSize(offs,iv,speci,tile_box));
pindxj = (int)floor(amrex::Random(engine)*getBinSize(offs,iv,specj,tile_box));
pindxi = specLists[speci][pindxi];
pindxj = specLists[specj][pindxj];
//pindxi = m_cell_vectors[speci][grid_id][imap][pindxi];
//pindxj = m_cell_vectors[specj][grid_id][imap][pindxj];

//if(i==1){Print() << "spec " << speci << " part " << pindxi << " of " << m_cell_vectors[speci][grid_id][imap].size() << endl;}
//if(i==1){Print() << "spec " << specj << " part " << pindxj << " of " << m_cell_vectors[specj][grid_id][imap].size() << endl;}
Expand All @@ -284,16 +274,15 @@ void FhdParticleContainer::CollideParticles(Real dt)
vj[1] = partj.rdata(FHD_realData::vely);
vj[2] = partj.rdata(FHD_realData::velz);


//if(i==1){Print() << "vel1 " << vi[0] << " vel2 " << vj[0] << endl;}

vij[0] = vi[0]-vj[0]; vij[1] = vi[1]-vj[1]; vij[2] = vi[2]-vj[2];
vrmag = sqrt(pow(vij[0],2)+pow(vij[1],2)+pow(vij[2],2));
if(vrmag>vrmax) {vrmax = vrmag; arrvrmax(i,j,k,specij) = 1.1*vrmax;}

Real eijmag = 0;
Real theta = 2.0*pi_usr*amrex::Random();
Real phi = std::acos(1.0-2.0*amrex::Random());
Real theta = 2.0*M_PI*amrex::Random(engine);
Real phi = std::acos(1.0-2.0*amrex::Random(engine));
eij[0] = std::sin(phi)*std::cos(theta);
eij[1] = std::sin(phi)*std::sin(theta);
eij[2] = std::cos(phi);
Expand All @@ -305,7 +294,7 @@ void FhdParticleContainer::CollideParticles(Real dt)
}

vreijmag = vij[0]*eij[0]+vij[1]*eij[1]+vij[2]*eij[2];
if(vrmag>vrmax*amrex::Random())
if(vrmag>vrmax*amrex::Random(engine))
{
vreijmag = vreijmag*2.0/massij;
vreij[0] = vreijmag*eij[0];
Expand All @@ -321,17 +310,11 @@ void FhdParticleContainer::CollideParticles(Real dt)
totalCol++;
}
}
for(int kk=0;kk<4;kk++)
{
if(i==1)
{
//Print() << "post ij2 " << kk << ": " << NSel[kk] << endl;
}
}
//Print() << "total cols: " << totalCol << "\n";
}
}
}

// }
// }
// }
});
}
}

Expand Down
20 changes: 15 additions & 5 deletions exec/DSMC/GNUmakefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ USE_OMP = FALSE
USE_CUDA = FALSE
COMP = gnu
DIM = 3
MAX_SPEC = 8

TINY_PROFILE = FALSE
USE_PARTICLES = TRUE
Expand Down Expand Up @@ -52,9 +53,9 @@ include $(AMREX_HOME)/Src/Boundary/Make.package
include $(AMREX_HOME)/Src/AmrCore/Make.package
include $(AMREX_HOME)/Src/Particle/Make.package

include $(AMREX_HOME)/Src/Extern/SWFFT/Make.package
INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT
VPATH_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT
#include $(AMREX_HOME)/Src/Extern/SWFFT/Make.package
#INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT
#VPATH_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT

#ifeq ($(NERSC_HOST),) # check if NERSC_HOST is empty => hack to test if on cori
# LIBRARIES += -L$(FFTW_DIR) -lfftw3
Expand All @@ -68,13 +69,22 @@ VPATH_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT
# $(info NERSC_HOST not empty => relying on compiler wrapper for fftw3)
#endif

LIBRARIES += -L$(FFTW_DIR) -lfftw3_mpi -lfftw3_omp -lfftw3 -lgomp
#LIBRARIES += -L$(FFTW_DIR) -lfftw3_mpi -lfftw3_omp -lfftw3 -lgomp
#LIBRARIES += -L$(FFTW_DIR) -lfftw3

include $(AMREX_HOME)/Tools/GNUMake/Make.rules

#DEFINES += -DDSMC=$(DSMC)

ifeq ($(findstring cgpu, $(HOST)), cgpu)
CXXFLAGS += $(FFTW)
CXXFLAGS += $(FFTW)
endif

ifeq ($(USE_CUDA),TRUE)
LIBRARIES += -lcufft
else
LIBRARIES += -L$(FFTW_DIR) -lfftw3_mpi -lfftw3
endif

MAXSPECIES := $(strip $(MAX_SPEC))
DEFINES += -DMAX_SPECIES=$(MAXSPECIES)
1 change: 1 addition & 0 deletions exec/DSMC/delPlots
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,3 +3,4 @@ rm -r pltpart*
rm -r pltprim*
rm -r pltcu*
rm -r spatialCross1D*
rm -r plt_SF_prim_*
Loading

0 comments on commit 7bebd68

Please sign in to comment.