capability for adding helix

Tests/Particles/Mapped/InitHelix.cpp

@@ -0,0 +1,55 @@
+#include <AMReX.H>
+#include <AMReX_MultiFab.H>
+
+using namespace amrex;
+
+void
+InitHelix (amrex::MultiFab& a_xyz_loc, amrex::Geometry& geom)
+{
+ const Real tpi = 2.0* amrex::Math::pi<Real>();
+
+ auto domain = geom.Domain();
+ auto probhi = geom.ProbHiArray();
+ auto problo = geom.ProbLoArray();
+
+ Real ilen = static_cast<Real>(domain.length(0));
+ Real jlen = static_cast<Real>(domain.length(1));
+ Real klen = static_cast<Real>(domain.length(2));
+
+ // Center of the toroidal quadrilateral
+ Real cx = 0.5 * (problo[0]+probhi[0]);
+ Real cy = 0.5 * (problo[1]+probhi[1]);
+ Real cz = 0.5 * (problo[2]+probhi[2]);
+
+ // This is "delta eta" which is in the j direction
+ Real d_eta = 1. / jlen;
+
+ // This is "delta xi" which is in the i direction
+ Real d_xi = 1. / ilen;
+
+ // This is "delta zeta" which is in the k direction
+ Real d_zeta = 1. / klen;
+ // loc_arr is nodal so no offset
+ for(MFIter mfi(a_xyz_loc); mfi.isValid(); ++mfi)
+ {
+ const Box& gtbx = mfi.growntilebox();
+ amrex::Print() << " grown tile box : " << gtbx << "\n";
+ auto loc_arr = a_xyz_loc.array(mfi);
+
+ ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ {
+ Real xi = (static_cast<int>(i)) * d_xi;
+ Real eta = (static_cast<int>(j)) * d_eta;
+ Real zeta = (static_cast<int>(k)) * d_zeta;
+
+ loc_arr(i,j,k,0) = cx + (0.1 * xi + 0.1) * cos(tpi * eta);
+ loc_arr(i,j,k,1) = cy + (0.1 * xi + 0.1) * sin(tpi * eta);
+ //loc_arr(i,j,k,2) = cz + ( xi * 0.1 - 0.1 ) + zeta * 0.1 * (2. + 1. * xi); 3D torus
+ loc_arr(i,j,k,2) = zeta*0.02 + 0.02*tpi*eta;// ( xi * 0.1 - 0.1 ) + zeta * 0.1 * (2. + 1. * xi) * tpi*eta;
+ amrex::Real rad = std::sqrt((loc_arr(i,j,k,0)-0.5)*(loc_arr(i,j,k,0)-0.5) + (loc_arr(i,j,k,1)-0.5)*(loc_arr(i,j,k,1)-0.5)) ;
+ });
+
+ }
+ a_xyz_loc.FillBoundary(geom.periodicity());
+
+}

Tests/Particles/Mapped/InitTorus.cpp

@@ -44,8 +44,7 @@ InitTorus (amrex::MultiFab& a_xyz_loc, amrex::Geometry& geom)
 
  loc_arr(i,j,k,0) = cx + (0.1 * xi + 0.1) * cos(tpi * eta);
  loc_arr(i,j,k,1) = cy + (0.1 * xi + 0.1) * sin(tpi * eta);
- //loc_arr(i,j,k,2) = cz + ( xi * 0.1 - 0.1 ) + zeta * 0.1 * (2. + 1. * xi); 3D torus
- loc_arr(i,j,k,2) = cz + zeta*0.2 - 0.1 + 0.3;
+ loc_arr(i,j,k,2) = cz + ( xi * 0.1 - 0.1 ) + zeta * 0.1 * (2. + 1. * xi); //3D torus
  amrex::Real rad = std::sqrt((loc_arr(i,j,k,0)-0.5)*(loc_arr(i,j,k,0)-0.5) + (loc_arr(i,j,k,1)-0.5)*(loc_arr(i,j,k,1)-0.5)) ;
  });
 

Tests/Particles/Mapped/InitUND.cpp

@@ -4,7 +4,7 @@
 using namespace amrex;
 
 enum struct ProbType {
- Torus, Annulus, Stretched, Hill
+ Helix, Torus, Annulus, Stretched, Hill
 };
 
 void
@@ -18,11 +18,50 @@ InitUND_map (MultiFab& und, const MultiFab& a_xyz_loc, Geometry& geom, int flow_
  // Center of the annulus
  Real cx = 0.5 * (problo[0]+probhi[0]);
  Real cy = 0.5 * (problo[1]+probhi[1]);
-
+ const Real tpi = 2.0* amrex::Math::pi<Real>();
  //
  // ANNULUS
  //
- if (prob_type == ProbType::Torus) {
+ if (prob_type == ProbType::Helix) {
+ // 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));
+ theta = std::atan2((y-cy),(x-cx));
+ und_arr(i,j,k,0) = -tpi* rad * std::sin(theta); // rad*sin(theta);
+ und_arr(i,j,k,1) = tpi * rad * std::cos(theta); // -1.*rad*cos(theta);
+
+#if (AMREX_SPACEDIM == 3)
+ // Real z = loc_arr(i,j,k,2);
+ und_arr(i,j,k,2) = 0.02*tpi;
+#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::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)

Tests/Particles/Mapped/Make.package

@@ -1,6 +1,7 @@
 CEXE_sources += main.cpp
 
 CEXE_sources += InitTorus.cpp
+CEXE_sources += InitHelix.cpp
 CEXE_sources += InitAnnulus.cpp
 CEXE_sources += InitHill.cpp
 CEXE_sources += InitStretched.cpp

Tests/Particles/Mapped/MappedPC.cpp

@@ -43,7 +43,7 @@ InitParticles (MultiFab& a_xyz_loc)
 #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) { // FOR EVERYTHING BUT TORUS
- if ( ( (iv[0] > 1) && (iv[0] < 30) ) && (iv[1] == 0 ) && (iv[2]>1 && iv[2] < 10 )) { // FOR TORUS
+ if ( ( (iv[0] > 4) && (iv[0] < 28) ) && (iv[1] == 1 ) && (iv[2]>1 && iv[2] < 50 )) { // FOR TORUS
 #endif
  int i = iv[0];
  int j = iv[1];

Tests/Particles/Mapped/inputs

@@ -2,11 +2,11 @@ size = (32,32,64 )
 max_grid_size = 128
 is_periodic = (1,1,1)
 num_ppc = 1
-nsteps = 750
+nsteps = 1
 nlevs = 1
 
 grid_type = mapped
 vel_type = nd
-prob_type = torus
+prob_type = helix
 
 do_plotfile = 1