30 PSRead(amrex::Vector<amrex::Geometry>& _geom) :
IC(_geom) {}
44 int ncomp = a_phi[lev]->nComp();
46 for (amrex::MFIter mfi(*a_phi[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
49 amrex::IndexType type = a_phi[lev]->ixType();
50 if (type == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
51 else if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
54 amrex::Array4<Set::Scalar>
const& phi = a_phi[lev]->array(mfi);
55 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k)
60 for (
unsigned int n = 0; n <
X.size(); n++)
63 Set::Scalar d1 = std::numeric_limits<Set::Scalar>::infinity();
65 d = (x -
X[n]).lpNorm<2>();
67 if (geom[0].isPeriodic(0))
69 d1 = std::min((x -
X[n] + size(0) * Set::Vector::Unit(0)).lpNorm<2>(),
70 (x -
X[n] - size(0) * Set::Vector::Unit(0)).lpNorm<2>());
72 #if AMREX_SPACEDIM > 1
73 if (geom[0].isPeriodic(1))
75 d1 = std::min((x -
X[n] + size(1) * Set::Vector::Unit(1)).lpNorm<2>(),
76 (x -
X[n] - size(1) * Set::Vector::Unit(1)).lpNorm<2>());
79 #if AMREX_SPACEDIM > 2
80 if (geom[0].isPeriodic(2))
82 d1 = std::min((x -
X[n] + size(2) * Set::Vector::Unit(2)).lpNorm<2>(),
83 (x -
X[n] - size(2) * Set::Vector::Unit(2)).lpNorm<2>());
89 if (d <= (R[n] + eps))
92 Set::Scalar m = 0.5 * (1 + erf((-d + R[n]) / (eps)));
93 min_grain_id = min_grain_id + m * (1. - min_grain_id);
97 phi(i, j, k, 0) = min_grain_id;
98 if (ncomp > 1) phi(i, j, k, 1) = 1.0 - min_grain_id;
115 if (datafile.is_open())
120 while (getline(datafile, line))
122 std::istringstream in(line);
124 std::string strx, stry, strz, strR;
125 in >> strx >> stry >> strz >> strR;
127 value.
R.push_back(value.
mult*std::stod(strR));
140 std::vector<Set::Vector>
X;
141 std::vector<Set::Scalar>
R;