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 : #include "Util/Util.H"
12 : using namespace std;
13 : #include <iostream>
14 : #include <vector>
15 : #include <algorithm>
16 : #include <iterator>
17 : #include "mpi.h"
18 : #include "IO/ParmParse.H"
19 :
20 : namespace IC
21 : {
22 : class PSRead : public IC<Set::Scalar>
23 : {
24 : public:
25 : static constexpr const char* name = "psread";
26 :
27 : enum Type
28 : {
29 : Partition,
30 : Values
31 : };
32 :
33 : PSRead(amrex::Vector<amrex::Geometry>& _geom) : IC(_geom) {}
34 :
35 0 : PSRead(amrex::Vector<amrex::Geometry>& _geom, IO::ParmParse& pp, std::string name) : IC(_geom)
36 : {
37 0 : pp.queryclass(name, *this);
38 0 : }
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 0 : }
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 : // Diffuseness of the sphere boundary
111 0 : pp.query_default("eps", value.eps, "0.0", Unit::Length());
112 :
113 0 : pp.forbid("filename", "use file.name instead");
114 : // Location of .xyzr file
115 0 : pp.query_file("file.name", filename);
116 :
117 : // Verbosity (used in parser only)
118 0 : pp.query_default("verbose", verbose, 0);
119 :
120 : // Unitless coordinate multiplier
121 0 : pp.query_default("mult",value.mult, "1.0", Unit::Less());
122 :
123 : // Whether to invert (1 outside instead of inside spheres)
124 0 : pp.query_default("invert",value.invert,false);
125 :
126 : // X-offset
127 0 : pp.queryarr_default("x0",value.x0,"0 0 0", Unit::Length());
128 :
129 0 : if (pp.contains("x0"))
130 0 : pp_queryarr("x0",value.x0); // Coordinate offset
131 :
132 0 : Unit unit = Unit::Length();
133 : // Units of length in the file
134 0 : pp.queryunit("file.unit",unit);
135 0 : Util::AssertException( INFO,TEST(unit.isType(Unit::Length())),
136 0 : "Unit must be of type length but got unit ",unit.normalized_unitstring());
137 :
138 0 : std::ifstream datafile(filename);
139 0 : std::string line;
140 0 : if (datafile.is_open())
141 : {
142 0 : value.X.clear();
143 0 : value.R.clear();
144 :
145 0 : while (getline(datafile, line))
146 : {
147 0 : std::istringstream in(line);
148 :
149 0 : std::string strx, stry, strz, strR;
150 0 : in >> strx >> stry >> strz >> strR;
151 :
152 0 : Set::Scalar x = (std::stod(strx) * unit).normalized_value();
153 0 : Set::Scalar y = (std::stod(stry) * unit).normalized_value();
154 : #if AMREX_SPACEDIM > 2
155 0 : Set::Scalar z = (std::stod(strz) * unit).normalized_value();
156 : #endif
157 0 : Set::Scalar r = (std::stod(strR) * unit).normalized_value();
158 :
159 0 : Set::Vector X(AMREX_D_DECL(x,y,z));
160 0 : X = value.x0 + value.mult*X;
161 :
162 0 : value.X.push_back(X);
163 0 : value.R.push_back(r);
164 0 : if (verbose > 0)
165 0 : Util::Message(INFO, "x=", value.X.back().transpose(), " r=", value.R.back());
166 0 : }
167 0 : datafile.close();
168 : }
169 : else
170 : {
171 0 : Util::Abort(INFO, "Unable to open file ", filename);
172 : }
173 0 : }
174 :
175 : private:
176 : std::vector<Set::Vector> X;
177 : std::vector<Set::Scalar> R;
178 : Set::Scalar eps;
179 : Set::Scalar mult = 1.0;
180 : Set::Vector x0 = Set::Vector::Zero();
181 : bool invert = false;
182 : };
183 : }
184 : #endif
|
