Line data    Source code

       1             : //
       2             : // Fill a domain (region where field=0) with packed spheres (regions where field=1).
       3             : // Sphere locations and radii are determined from an xyzr file.
       4             : //
       5             : 
       6             : #ifndef IC_PSREAD_H_
       7             : #define IC_PSREAD_H_
       8             : 
       9             : #include "Set/Set.H"
      10             : #include "IC/IC.H"
      11             : using namespace std;
      12             : #include <iostream>
      13             : #include <vector>
      14             : #include <algorithm>
      15             : #include <iterator>
      16             : #include "mpi.h"
      17             : #include "IO/ParmParse.H"
      18             : 
      19             : namespace IC
      20             : {
      21             : class PSRead : public IC
      22             : {
      23             : public:
      24             :     static constexpr const char* name = "psread";
      25             : 
      26             :     enum Type
      27             :     {
      28             :         Partition,
      29             :         Values
      30             :     };
      31             : 
      32             :     PSRead(amrex::Vector<amrex::Geometry>& _geom) : IC(_geom) {}
      33             : 
      34           0 :     PSRead(amrex::Vector<amrex::Geometry>& _geom, IO::ParmParse& pp, std::string name) : IC(_geom)
      35             :     {
      36           0 :         pp_queryclass(name, *this);
      37           0 :     }
      38             : 
      39             :     void Define() {
      40             : 
      41             :     };
      42             : 
      43           0 :     void Add(const int& lev, Set::Field<Set::Scalar>& a_phi, Set::Scalar)
      44             :     {
      45           0 :         Set::Vector size = Set::Size(geom[lev]);
      46           0 :         int ncomp = a_phi[lev]->nComp();
      47             : 
      48           0 :         for (amrex::MFIter mfi(*a_phi[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
      49             :         {
      50           0 :             amrex::Box bx;
      51           0 :             amrex::IndexType type = a_phi[lev]->ixType();
      52           0 :             if (type == amrex::IndexType::TheCellType())      bx = mfi.growntilebox();
      53           0 :             else if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
      54           0 :             else Util::Abort(INFO, "Unkonwn index type");
      55             : 
      56           0 :             amrex::Array4<Set::Scalar> const& phi = a_phi[lev]->array(mfi);
      57           0 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
      58             :             {
      59           0 :                 Set::Vector x = Set::Position(i, j, k, geom[lev], type);
      60           0 :                 Set::Scalar min_grain_id = 0;
      61             : 
      62           0 :                 for (unsigned int n = 0; n < X.size(); n++)
      63             :                 {
      64             :                     Set::Scalar d;
      65           0 :                     Set::Scalar d1 = std::numeric_limits<Set::Scalar>::infinity();
      66             : 
      67           0 :                     d = (x - X[n]).lpNorm<2>();
      68             : 
      69           0 :                     if (geom[0].isPeriodic(0))
      70             :                     {
      71           0 :                         d1 = std::min((x - X[n] + size(0) * Set::Vector::Unit(0)).lpNorm<2>(),
      72           0 :                             (x - X[n] - size(0) * Set::Vector::Unit(0)).lpNorm<2>());
      73             :                     }
      74             : #if AMREX_SPACEDIM > 1
      75           0 :                     if (geom[0].isPeriodic(1))
      76             :                     {
      77           0 :                         d1 = std::min((x - X[n] + size(1) * Set::Vector::Unit(1)).lpNorm<2>(),
      78           0 :                             (x - X[n] - size(1) * Set::Vector::Unit(1)).lpNorm<2>());
      79             :                     }
      80             : #endif
      81             : #if AMREX_SPACEDIM > 2
      82           0 :                     if (geom[0].isPeriodic(2))
      83             :                     {
      84           0 :                         d1 = std::min((x - X[n] + size(2) * Set::Vector::Unit(2)).lpNorm<2>(),
      85           0 :                             (x - X[n] - size(2) * Set::Vector::Unit(2)).lpNorm<2>());
      86             :                     }
      87             : #endif
      88             : 
      89           0 :                     d = std::min(d, d1);
      90             : 
      91           0 :                     if (d <= (R[n] + eps))
      92             : 
      93             :                     {
      94           0 :                         Set::Scalar m = 0.5 * (1 + erf((-d + R[n]) / (eps)));
      95           0 :                         min_grain_id = min_grain_id + m * (1. - min_grain_id);
      96             :                     }
      97             :                 }
      98             : 
      99           0 :                 phi(i, j, k, 0) = min_grain_id;
     100           0 :                 if (invert) phi(i,j,k,0) = 1.0 - min_grain_id;
     101           0 :                 if (ncomp > 1) phi(i, j, k, 1) = 1.0 - min_grain_id;
     102           0 :             });
     103             :         }
     104           0 :     }
     105             : 
     106           0 :     static void Parse(PSRead& value, IO::ParmParse& pp)
     107             :     {
     108           0 :         std::string filename;
     109           0 :         int verbose = 0;
     110           0 :         pp_query("eps", value.eps); // Diffuseness of the sphere boundary
     111           0 :         pp_query_file("filename", filename); // Location of .xyzr file
     112           0 :         pp_query("verbose", verbose); // Verbosity (used in parser only)
     113           0 :         pp_query("mult",value.mult); // Coordinate multiplier
     114           0 :         pp.query("invert",value.invert); // Coordinate multiplier
     115           0 :         if (pp.contains("x0"))
     116           0 :             pp_queryarr("x0",value.x0); // Coordinate offset
     117           0 :         std::ifstream datafile(filename);
     118           0 :         std::string line;
     119           0 :         if (datafile.is_open())
     120             :         {
     121           0 :             value.X.clear();
     122           0 :             value.R.clear();
     123             : 
     124           0 :             while (getline(datafile, line))
     125             :             {
     126           0 :                 std::istringstream in(line);
     127             : 
     128           0 :                 std::string strx, stry, strz, strR;
     129           0 :                 in >> strx >> stry >> strz >> strR;
     130           0 :                 value.X.push_back(value.x0 + value.mult*Set::Vector(AMREX_D_DECL(std::stod(strx), std::stod(stry), std::stod(strz))));
     131           0 :                 value.R.push_back(value.mult*std::stod(strR));
     132           0 :                 if (verbose > 0)
     133           0 :                     Util::Message(INFO, "x=", value.X.back().transpose(), " r=", value.R.back());
     134             :             }
     135           0 :             datafile.close();
     136             :         }
     137             :         else
     138             :         {
     139           0 :             Util::Abort(INFO, "Unable to open file ", filename);
     140             :         }
     141           0 :     }
     142             : 
     143             : private:
     144             :     std::vector<Set::Vector> X;
     145             :     std::vector<Set::Scalar> R;
     146             :     Set::Scalar eps;
     147             :     Set::Scalar mult = 1.0;
     148             :     Set::Vector x0 = Set::Vector::Zero();
     149             :     bool invert = false;
     150             : };
     151             : }
     152             : #endif