Skip to content

Commit

Permalink
3d particle mapping
Browse files Browse the repository at this point in the history
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Feb 3, 2024
1 parent 7393481 commit d413ba5
Show file tree
Hide file tree
Showing 5 changed files with 138 additions and 42 deletions.
43 changes: 41 additions & 2 deletions Tests/Particles/Mapped/InitUND.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
using namespace amrex;

enum struct ProbType {
Annulus, Stretched, Hill
Torus, Annulus, Stretched, Hill
};

void
Expand All @@ -22,7 +22,46 @@ InitUND_map (MultiFab& und, const MultiFab& a_xyz_loc, Geometry& geom, int flow_
//
// ANNULUS
//
if (prob_type == ProbType::Annulus) {
if (prob_type == ProbType::Torus) {
// We only do this problem with fully mapped coordinates
AMREX_ALWAYS_ASSERT(a_xyz_loc.nComp() == AMREX_SPACEDIM);
for (MFIter mfi(und); mfi.isValid(); ++mfi)
{
const Box& tile_box = mfi.growntilebox();
auto loc_arr = a_xyz_loc.const_array(mfi);
auto und_arr = und.array(mfi);

ParallelFor( tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// Physical location of cell center
Real x = loc_arr(i,j,k,0);
Real y = loc_arr(i,j,k,1);
Real theta;
//if (x == cx) {
// theta = 3.14/2.;
//} else {
// theta = atan((y-cy)/(x-cx));
//}
//Real rad = sqrt( (x-cx)*( x-cx) + (y-cy)*(y-cy));

und_arr(i,j,k,0) = y - cy; // rad*sin(theta);
und_arr(i,j,k,1) = cx - x; // -1.*rad*cos(theta);

#if (AMREX_SPACEDIM == 3)
// Real z = loc_arr(i,j,k,2);
und_arr(i,j,k,2) = 0.02;
#endif
if (i == 20 && j==5 && k==5) amrex::Print() << "UND AT " << IntVect(AMREX_D_DECL(i,j,k)) << " "
<< RealVect(AMREX_D_DECL(und_arr(i,j,k,0),und_arr(i,j,k,1),und_arr(i,j,k,2)))
<< std::endl;
if (i == 20 && j==23 && k==5) amrex::Print() << "UND AT " << IntVect(AMREX_D_DECL(i,j,k)) << " "
<< RealVect(AMREX_D_DECL(und_arr(i,j,k,0),und_arr(i,j,k,1),und_arr(i,j,k,2)))
<< std::endl;


});
}
} else if (prob_type == ProbType::Annulus) {

// We only do this problem with fully mapped coordinates
AMREX_ALWAYS_ASSERT(a_xyz_loc.nComp() == AMREX_SPACEDIM);
Expand Down
70 changes: 56 additions & 14 deletions Tests/Particles/Mapped/MappedPC.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,6 @@ public:
} else { neg++;
}
}
amrex::Print() << " pos " << pos << " neg " << neg << "\n";
inside = ( (pos > 0) && (neg > 0) ) ? false : true;
return inside;
}
Expand All @@ -113,7 +112,7 @@ public:
amrex::Array2D<int, 0, 3, 0, 5> inode{4, 2, 6, 0, 7, 1, 5, 3,
5, 0, 4, 1, 6, 3, 7, 2,
2, 1, 3, 0, 7, 4, 6, 5};
bool inside;
bool inside = false;
// Define ray 1 : particle to lower corner
amrex::Real r1x = x[0] - xp;
amrex::Real r1y = y[0] - yp;
Expand All @@ -126,15 +125,16 @@ public:
int numfaces = 6;
int numinside_counter = 0;
for (int i = 0; i < numfaces; ++i) {
amrex::Real nx = 0.5 * (y[inode(0,i)] - y[inode(1,i)]) * (z[inode(2,i)] - z[inode(3,i)])
- (z[inode(0,i)] - z[inode(1,i)]) * (y[inode(2,i)] - y[inode(3,i)]);
amrex::Real ny = 0.5 * (z[inode(0,i)] - z[inode(1,i)]) * (x[inode(2,i)] - x[inode(3,i)])
- (x[inode(0,i)] - x[inode(1,i)]) * (z[inode(2,i)] - z[inode(3,i)]);
amrex::Real nz = 0.5 * (x[inode(0,i)] - x[inode(1,i)]) * (y[inode(2,i)] - y[inode(3,i)])
- (y[inode(0,i)] - y[inode(1,i)]) * (x[inode(2,i)] - x[inode(3,i)]);
amrex::Real raydotn = (i%2 == 0) ? r1x * nx + r1y * ny + r1z * nz
: r2x * nx + r2y * ny + r2z * nz;
amrex::Real nx = 0.5 * ( ( (y[inode(0,i)] - y[inode(1,i)]) * (z[inode(2,i)] - z[inode(3,i)]) )
- ( (z[inode(0,i)] - z[inode(1,i)]) * (y[inode(2,i)] - y[inode(3,i)]) ) );
amrex::Real ny = 0.5 * ( ( (z[inode(0,i)] - z[inode(1,i)]) * (x[inode(2,i)] - x[inode(3,i)]) )
- ( (x[inode(0,i)] - x[inode(1,i)]) * (z[inode(2,i)] - z[inode(3,i)]) ) );
amrex::Real nz = 0.5 * ( ( (x[inode(0,i)] - x[inode(1,i)]) * (y[inode(2,i)] - y[inode(3,i)]) )
- ( (y[inode(0,i)] - y[inode(1,i)]) * (x[inode(2,i)] - x[inode(3,i)]) ) );
amrex::Real raydotn = (i%2 == 0) ? (r1x * nx + r1y * ny + r1z * nz)
: (r2x * nx + r2y * ny + r2z * nz);
if (raydotn < 0) {
amrex::Print() << " iface : " << i << " outside \n";
inside = false;
return inside;
}
Expand Down Expand Up @@ -195,6 +195,7 @@ public:
}}
}

if (!inside) amrex::Abort("particle has not been mapped in 2D!\n");
#elif (AMREX_SPACEDIM == 3)
//amrex::IntVect iv( int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
// int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
Expand All @@ -218,9 +219,9 @@ public:
//}

amrex::IntVect iv( p.idata(0), p.idata(1), p.idata(2) );
amrex::Real x[8] = {0, 1, 0, 1, 0, 1, 0, 1};
amrex::Real y[8] = {0, 0, 1, 1, 0, 0, 1, 1};
amrex::Real z[8] = {0, 0, 0, 0, 1, 1, 1, 1};
amrex::Real x[8];
amrex::Real y[8];
amrex::Real z[8];

int cnt = 0;
for (int k = 0; k < 2; ++k) {
Expand Down Expand Up @@ -256,11 +257,52 @@ public:
p.idata(0) = p.idata(0) + lo_i;
p.idata(1) = p.idata(1) + lo_j;
p.idata(2) = p.idata(2) + lo_k;
if (p.idata(0) == -1) p.idata(0) = 31;
if (p.idata(1) == -1) p.idata(1) = 31;
if (p.idata(0) == 32) p.idata(0) = 0;
if (p.idata(1) == 32) p.idata(1) = 0;
return;
}
}}}
}
if (!inside) amrex::Abort("particle has not been mapped!\n");
if (!inside) {
// // looping over neighboring cells
// for (int lo_k = -2; lo_k <=2; ++lo_k){
// for (int lo_j = -2; lo_j <=2; ++lo_j){
// for (int lo_i = -2; lo_i <=2; ++lo_i){
// if ( (lo_i == 0) && (lo_j == 0) && lo_k == 0 ) continue;
// amrex::IntVect ivn(p.idata(0) + lo_i, p.idata(1) + lo_j, p.idata(2) + lo_k);
// cnt = 0;
// for (int k = 0; k < 2; ++k) {
// for (int j = 0; j < 2; ++j) {
// for (int i = 0; i < 2; ++i) {
// x[cnt] = loc_arr(ivn[0] + i, ivn[1] + j, ivn[2] + k, 0);
// y[cnt] = loc_arr(ivn[0] + i, ivn[1] + j, ivn[2] + k, 1);
// z[cnt] = loc_arr(ivn[0] + i, ivn[1] + j, ivn[2] + k, 2);
// cnt++;
// }}}
// amrex::Print() << " ngbr cnt : " << lo_i << " " << lo_j << " " << lo_k << " " << ivn<< " " << "\n";
// for (int ii = 0; ii < 8; ++ii ) {amrex::Print() << " " << x[ii] << ", ";}
// amrex::Print() << "\n";
// for (int ii = 0; ii < 8; ++ii ) {amrex::Print() << " " << y[ii] << ", ";}
// amrex::Print() << "\n";
// for (int ii = 0; ii < 8; ++ii ) {amrex::Print() << " " << z[ii] << ", ";}
// amrex::Print() << "\n";
// inside = PointInHexahedron (x, y, z, p.pos(0), p.pos(1), p.pos(2));
// if (inside) {
// p.idata(0) = p.idata(0) + lo_i;
// p.idata(1) = p.idata(1) + lo_j;
// p.idata(2) = p.idata(2) + lo_k;
// if (p.idata(0) == -1) p.idata(0) = 31;
// if (p.idata(1) == -1) p.idata(1) = 31;
// if (p.idata(0) == 32) p.idata(0) = 0;
// if (p.idata(1) == 32) p.idata(1) = 0;
// return;
// }
// }}}
// amrex::Print() << " pos : " << p.pos(0) << " " << p.pos(1) << " " << p.pos(2) << " iv " << iv << "\n";
if (!inside) amrex::Abort("particle has not been mapped!\n");
}
#endif
}

Expand Down
49 changes: 33 additions & 16 deletions Tests/Particles/Mapped/MappedPC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,23 +42,27 @@ InitParticles (MultiFab& a_xyz_loc)
if (iv[0] == 0) { // FOR ANNULUS
#elif (AMREX_SPACEDIM == 3)
int k = iv[2];
if (iv[0] == 0 && iv[1] == 0 && iv[2] >= dom_lo.z && iv[2] <= dom_hi.z/2) {
//if (iv[0] == 0 && iv[1] == 0 && iv[2] >= dom_lo.z && iv[2] <= dom_hi.z/2) { // FOR EVERYTHING BUT TORUS
if ( ( (iv[0] > 1) && (iv[0] < 30) ) && (iv[1] == 0 ) && (iv[2]>1 && iv[2] < 10 )) { // FOR TORUS
#endif
int i = iv[0];
int j = iv[1];

// This is the physical location of the center of the cell
Real x = 0.25*( loc_arr(i ,j,k,0) + loc_arr(i ,j+1,k,0)
+loc_arr(i+1,j,k,0) + loc_arr(i+1,j+1,k,0));
Real y = 0.25*( loc_arr(i ,j,k,1) + loc_arr(i ,j+1,k,1)
+loc_arr(i+1,j,k,1) + loc_arr(i+1,j+1,k,1));
Real x = 0.125*( loc_arr(i ,j ,k ,0) + loc_arr(i+1,j ,k , 0)
+ loc_arr(i ,j+1,k ,0) + loc_arr(i+1,j+1,k , 0)
+ loc_arr(i ,j ,k+1,0) + loc_arr(i+1,j ,k+1, 0)
+ loc_arr(i ,j+1,k+1,0) + loc_arr(i+1,j+1,k+1, 0));
Real y = 0.125*( loc_arr(i ,j ,k ,1) + loc_arr(i+1,j ,k , 1)
+ loc_arr(i ,j+1,k ,1) + loc_arr(i+1,j+1,k , 1)
+ loc_arr(i ,j ,k+1,1) + loc_arr(i+1,j ,k+1, 1)
+ loc_arr(i ,j+1,k+1,1) + loc_arr(i+1,j+1,k+1, 1));

ParticleType p;
p.id() = ParticleType::NextID();
p.cpu() = ParallelDescriptor::MyProc();
p.pos(0) = x;
p.pos(1) = y;

p.rdata(MappedRealIdx::vx) =Real(0.0);
p.rdata(MappedRealIdx::vy) = Real(0.0);

Expand All @@ -68,15 +72,25 @@ InitParticles (MultiFab& a_xyz_loc)
#if (AMREX_SPACEDIM == 2)
amrex::Print() << "Particle at (x,y) OF " << iv << " " << x << " " << y << std::endl;
#elif (AMREX_SPACEDIM == 3)
Real z = Real(0.125)*(loc_arr(i ,j,k ,2) + loc_arr(i ,j+1,k ,2) +
loc_arr(i+1,j,k ,2) + loc_arr(i+1,j+1,k ,2) +
loc_arr(i ,j,k+1,2) + loc_arr(i ,j+1,k+1,2) +
loc_arr(i+1,j,k+1,2) + loc_arr(i+1,j+1,k+1,2));
Real z = 0.125*( loc_arr(i ,j ,k ,2) + loc_arr(i+1,j ,k , 2)
+ loc_arr(i ,j+1,k ,2) + loc_arr(i+1,j+1,k , 2)
+ loc_arr(i ,j ,k+1,2) + loc_arr(i+1,j ,k+1, 2)
+ loc_arr(i ,j+1,k+1,2) + loc_arr(i+1,j+1,k+1, 2));
p.pos(2) = z;
p.rdata(MappedRealIdx::vz) = Real(0.);
p.idata(MappedIntIdx::k) = iv[2]; // particles carry their k-index

amrex::Print() << "Particle at (x,y,z) OF " << iv << " " << x << " " << y << " " << z << std::endl;

// amrex::Print() << " i, j, k " << loc_arr(i ,j ,k ,0) << " " << loc_arr(i ,j ,k ,1) << " " << loc_arr(i ,j ,k ,2) << "\n";
// amrex::Print() << " i+1, j, k " << loc_arr(i+1 ,j ,k ,0) << " " << loc_arr(i+1 ,j ,k ,1) << " " << loc_arr(i+1 ,j ,k ,2) << "\n";
// amrex::Print() << " i+1, j+1, k " << loc_arr(i+1 ,j+1 ,k ,0) << " " << loc_arr(i+1 ,j+1 ,k ,1) << " " << loc_arr(i+1 ,j+1 ,k ,2) << "\n";
// amrex::Print() << " i, j+1, k " << loc_arr(i ,j+1 ,k ,0) << " " << loc_arr(i ,j+1 ,k ,1) << " " << loc_arr(i ,j+1 ,k ,2) << "\n";
// amrex::Print() << " i, j, k+1 " << loc_arr(i ,j ,k+1 ,0) << " " << loc_arr(i ,j ,k+1 ,1) << " " << loc_arr(i ,j ,k+1 ,2) << "\n";
// amrex::Print() << " i+1, j, k+1 " << loc_arr(i+1 ,j ,k+1 ,0) << " " << loc_arr(i+1 ,j ,k+1 ,1) << " " << loc_arr(i+1 ,j ,k+1 ,2) << "\n";
// amrex::Print() << " i+1, j+1, k+1 " << loc_arr(i+1 ,j+1 ,k+1 ,0) << " " << loc_arr(i+1 ,j+1 ,k+1 ,1) << " " << loc_arr(i+1 ,j+1 ,k+1 ,2) << "\n";
// amrex::Print() << " i, j+1, k+1 " << loc_arr(i ,j+1 ,k+1 ,0) << " " << loc_arr(i ,j+1 ,k+1 ,1) << " " << loc_arr(i ,j+1 ,k+1 ,2) << "\n";
// amrex::AllPrintToFile("Initpartpos.txt") << p.pos(0) << " " << p.pos(1) << " " << p.pos(2)<< " \n";
#endif

host_particles.push_back(p);
Expand Down Expand Up @@ -234,33 +248,36 @@ MappedPC::AdvectWithUND (MultiFab& vel_nd, int lev, Real dt, const MultiFab& a_x
#if (AMREX_SPACEDIM == 2)
amrex::Print() << "FROM " << p.pos(0) << " " << p.pos(AMREX_SPACEDIM-1) << std::endl;
#elif (AMREX_SPACEDIM == 3)
amrex::Print() << "FROM " << p.pos(0) << " " << p.pos(1) << " " << p.pos(AMREX_SPACEDIM-1) << std::endl;
amrex::Real rc = sqrt((p.pos(0)-0.5)*(p.pos(0)-0.5) + (p.pos(1)-0.5)*(p.pos(1)-0.5));
amrex::Real theta = std::atan2( (p.pos(2) - 0.5) , (p.pos(1) - 0.5));
amrex::Print() << "FROM " << p.pos(0) << " " << p.pos(1) << " " << p.pos(AMREX_SPACEDIM-1) << " " << rc << " theta : " << theta << std::endl;
#endif
for (int dim=0; dim < AMREX_SPACEDIM; dim++)
{
p.rdata(dim) = p.pos(dim);
p.pos(dim) += static_cast<ParticleReal>(ParticleReal(0.5)*dt*v[dim]);
}
update_mapped_idata(p,plo,dxi,loc_arr);
amrex::AllPrintToFile("Frompartpos.txt") << p.pos(0) << " " << p.pos(1) << " " << p.pos(2)<< " \n";
// amrex::AllPrintToFile("Frompartpos.txt") << p.pos(0) << " " << p.pos(1) << " " << p.pos(2)<< " \n";
}
else
{
for (int dim=0; dim < AMREX_SPACEDIM; dim++)
{
p.pos(dim) = p.rdata(dim) + static_cast<ParticleReal>(dt*v[dim]);
p.rdata(dim) = v[dim];
amrex::Print() << dim << " pos : " << p.pos(dim) << " rdata : " << v[dim] << "\n";
}
#if (AMREX_SPACEDIM == 2)
amrex::Print() << "TO " << p.pos(0) << " " << p.pos(AMREX_SPACEDIM-1) <<
" WITH VEL " << v[0] << std::endl;
#elif (AMREX_SPACEDIM == 3)
amrex::Print() << "TO " << p.pos(0) << " " << p.pos(1) << " " << p.pos(AMREX_SPACEDIM-1) <<
" WITH VEL " << v[0] << std::endl;
amrex::Real rc = sqrt((p.pos(0)-0.5)*(p.pos(0)-0.5) + (p.pos(1)-0.5)*(p.pos(1)-0.5));
amrex::Real theta = std::atan2( (p.pos(2) - 0.5) , (p.pos(1) - 0.5));
amrex::Print() << "TO " << p.pos(0) << " " << p.pos(1) << " " << p.pos(AMREX_SPACEDIM-1) << " " << rc << " " << theta <<
" WITH VEL " << v[0] << " " << v[1] << " " << v[2] << std::endl;
#endif
update_mapped_idata(p,plo,dxi,loc_arr);
amrex::AllPrintToFile("partpos.txt") << p.pos(0) << " " << p.pos(1) << " " << p.pos(2)<< " \n";
// amrex::AllPrintToFile("partpos.txt") << p.pos(0) << " " << p.pos(1) << " " << p.pos(2)<< " \n";
}
});
} // ParIter
Expand Down
6 changes: 3 additions & 3 deletions Tests/Particles/Mapped/inputs
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
size = (32, 32, 32)
size = (32,32,64 )
max_grid_size = 128
is_periodic = (1,0)
is_periodic = (1,1,1)
num_ppc = 1
nsteps = 1
nsteps = 750
nlevs = 1

grid_type = mapped
Expand Down
12 changes: 5 additions & 7 deletions Tests/Particles/Mapped/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,12 +123,12 @@ void get_test_params(TestParams& params)

std::string prob_type_string;
pp.get("prob_type", prob_type_string);
AMREX_ALWAYS_ASSERT(prob_type_string == "donut" ||
AMREX_ALWAYS_ASSERT(prob_type_string == "torus" ||
prob_type_string == "annulus" ||
prob_type_string == "stretched" ||
prob_type_string == "unstretched" ||
prob_type_string == "hill");
if (prob_type_string == "donut" ) params.prob_type = ProbType::Torus;
if (prob_type_string == "torus" ) params.prob_type = ProbType::Torus;
if (prob_type_string == "annulus" ) params.prob_type = ProbType::Annulus;
if (prob_type_string == "unstretched") params.prob_type = ProbType::Unstretched;
if (prob_type_string == "stretched" ) params.prob_type = ProbType::Stretched;
Expand Down Expand Up @@ -178,7 +178,7 @@ void Test()

#if (AMREX_SPACEDIM == 3)
real_box.setLo(2, 0.0);
real_box.setHi(2, 1.0);
real_box.setHi(2, 5.0);
#endif

IntVect domain_lo(AMREX_D_DECL(0, 0, 0));
Expand Down Expand Up @@ -220,7 +220,7 @@ void Test()
MultiFab a_z_loc(ba_nd,dm[lev],1,1);

// This has AMREX_SPACEDIM components to define (x,y,z) as a function of (i,j,k)
MultiFab a_xyz_loc(ba_nd,dm[lev],AMREX_SPACEDIM,1);
MultiFab a_xyz_loc(ba_nd,dm[lev],AMREX_SPACEDIM,2);

// Annulus
if (params.prob_type == ProbType::Torus) {
Expand Down Expand Up @@ -278,7 +278,7 @@ void Test()
#endif

// Hard-wire the vertical velocity
Real vert_vel = 0.0;
Real vert_vel = 0.1;

if (params.vel_type == VelType::mac)
{
Expand Down Expand Up @@ -379,8 +379,6 @@ void Test()
amrex::Print() << "COMPUTING DT TO BE " << dt << " BASED ON MAX VEL " << max_vel << std::endl;
#else
auto dx = geom[0].CellSize();
amrex::Print() << dx[0] << " " << dx[1] << " " << dx[2] << "\n";
amrex::Print() << und.max(0,0,false) << "\n";
amrex::Real dt = 0.01;
amrex::Print() << "SETTING DT TO BE " << dt << std::endl;
#endif
Expand Down

0 comments on commit d413ba5

Please sign in to comment.