Line data Source code
1 : //
2 : // Initialize a sphere region in which the field \(\eta=1\).
3 : // Can be a 3D sphere (xyz) or oriented around any of the three
4 : // axes.
5 : //
6 : // :bdg-warning:`This functionality is replaced by IC::Expression.`
7 : //
8 :
9 : #ifndef IC_SPHERE_H_
10 : #define IC_SPHERE_H_
11 :
12 : #include <cmath>
13 :
14 : #include "IO/ParmParse.H"
15 : #include "IC/IC.H"
16 : #include "Unit/Unit.H"
17 :
18 : /// Initialize a spherical inclusion
19 : namespace IC
20 : {
21 : class Sphere: public IC<Set::Scalar>
22 : {
23 : public:
24 : static constexpr const char* name = "sphere";
25 :
26 : /// Sphere is also used for a cylindrical case, but switching from the default type to the desired orientation.
27 : /// For instance, XY refers to a cylinder with the circular faces in the XY plane.
28 : enum Type
29 : {
30 : XYZ,
31 : YZ,
32 : ZX,
33 : XY
34 : };
35 :
36 : Sphere(amrex::Vector<amrex::Geometry>& _geom): IC(_geom) {}
37 : Sphere(amrex::Vector<amrex::Geometry>& _geom, IO::ParmParse& pp): IC(_geom)
38 : {
39 : pp_queryclass(*this);
40 : }
41 1 : Sphere(amrex::Vector<amrex::Geometry>& _geom, IO::ParmParse& pp, std::string name): IC(_geom)
42 : {
43 1 : pp_queryclass(name, *this);
44 1 : }
45 2 : Sphere(amrex::Vector<amrex::Geometry>& _geom, Unit a_unit, IO::ParmParse& pp, std::string name):
46 2 : IC(_geom), unit(a_unit)
47 : {
48 2 : pp_queryclass(name, *this);
49 2 : }
50 : /// Constructor defining radius, center, dimension/orientation, and field values within and outside of the area.
51 : Sphere(amrex::Vector<amrex::Geometry>& _geom, Set::Scalar _radius, Set::Vector _center, Type _type = Type::XYZ, Set::Scalar _alpha_in = 1, Set::Scalar _alpha_out = 0)
52 : : IC(_geom)
53 : {
54 : Define(_radius, _center, _type, _alpha_in, _alpha_out);
55 : }
56 :
57 : void Define(Set::Scalar a_radius,
58 : Set::Vector a_center,
59 : Type a_type,
60 : Set::Scalar a_alpha_in,
61 : Set::Scalar a_alpha_out)
62 : {
63 : radius = a_radius;
64 : center = a_center;
65 : type = a_type;
66 : alpha_in = a_alpha_in;
67 : alpha_out = a_alpha_out;
68 : }
69 :
70 38 : void Add(const int& lev, Set::Field<Set::Scalar>& a_field, Set::Scalar)
71 : {
72 38 : bool cellcentered = (a_field[0]->boxArray().ixType() == amrex::IndexType(amrex::IntVect::TheCellVector()));
73 38 : int ncomp = a_field[0]->nComp();
74 :
75 3258 : for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
76 : {
77 3220 : amrex::Box bx = mfi.tilebox();
78 : // bx.grow(a_field[lev]->nGrow());
79 :
80 3220 : amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
81 3220 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
82 : {
83 :
84 : Set::Scalar AMREX_D_DECL(x, y, z);
85 207936 : if (cellcentered)
86 : {
87 207936 : AMREX_D_TERM(x = geom[lev].ProbLo()[0] + ((amrex::Real)(i)+0.5) * geom[lev].CellSize()[0];,
88 : y = geom[lev].ProbLo()[1] + ((amrex::Real)(j)+0.5) * geom[lev].CellSize()[1];,
89 : z = geom[lev].ProbLo()[2] + ((amrex::Real)(k)+0.5) * geom[lev].CellSize()[2];);
90 : }
91 : else
92 : {
93 0 : AMREX_D_TERM(x = geom[lev].ProbLo()[0] + (amrex::Real)(i)*geom[lev].CellSize()[0];,
94 : y = geom[lev].ProbLo()[1] + (amrex::Real)(j)*geom[lev].CellSize()[1];,
95 : z = geom[lev].ProbLo()[2] + (amrex::Real)(k)*geom[lev].CellSize()[2];);
96 : }
97 :
98 207936 : Set::Scalar rsq = NAN;
99 :
100 207936 : if (type == Type::XYZ)
101 : {
102 : // 3D Sphere
103 207936 : rsq =
104 207936 : AMREX_D_TERM((x - center(0)) * (x - center(0)),
105 : +(y - center(1)) * (y - center(1)),
106 : +(z - center(2)) * (z - center(2)));
107 : }
108 0 : else if (type == Type::XY)
109 : {
110 : // Cylinder along Z axis
111 0 : rsq =
112 0 : AMREX_D_TERM((x - center(0)) * (x - center(0)),
113 : +(y - center(1)) * (y - center(1)),
114 : );
115 : }
116 0 : else if (type == Type::YZ)
117 : {
118 : // Cylinder along X axis
119 0 : rsq =
120 0 : AMREX_D_TERM(,
121 : +(y - center(1)) * (y - center(1)),
122 : +(z - center(2)) * (z - center(2)));
123 : }
124 0 : else if (type == Type::ZX)
125 : {
126 : // Cylinder along Y axis
127 0 : rsq =
128 0 : AMREX_D_TERM((x - center(0)) * (x - center(0)),
129 : ,
130 : +(z - center(2)) * (z - center(2)));
131 : }
132 :
133 207936 : if (rsq < radius * radius)
134 : {
135 22968 : field(i, j, k, 0) = alpha_in;
136 45888 : if (ncomp > 1) field(i, j, k, 1) = 0.;
137 : }
138 : else
139 : {
140 184968 : field(i, j, k, 0) = alpha_out;
141 366336 : if (ncomp > 1) field(i, j, k, 1) = 1;
142 : }
143 207936 : });
144 38 : }
145 38 : };
146 : // This is a somewhat antiquated IC that will eventually be replaced
147 : // with the Expression IC.
148 3 : static void Parse(Sphere& value, IO::ParmParse& pp)
149 : {
150 : // Radius of the sphere
151 12 : pp_query_default("radius", value.radius,"1.0", Unit::Length());
152 : // Vector location of the sphere center
153 9 : pp_queryarr("center", value.center, Unit::Length());
154 : // Value of the field inside the sphere
155 12 : pp_query_default("inside", value.alpha_in,"1.0",value.unit);
156 : // Value of the field outside teh sphere
157 9 : pp_query_default("outside", value.alpha_out,"0.0", value.unit);
158 3 : std::string type;
159 : // Type - can be cylinder oriented along the x, y, z directions or full sphere.
160 12 : pp_query_validate("type", type, {"xyz","yz","zx","xy"});
161 3 : if (type == "yz") value.type = Type::YZ;
162 3 : if (type == "zx") value.type = Type::ZX;
163 3 : if (type == "xy") value.type = Type::XY;
164 3 : if (type == "xyz") value.type = Type::XYZ;
165 3 : }
166 :
167 : private:
168 : Set::Vector center = Set::Vector::Zero();
169 : Set::Scalar radius = NAN;
170 : Set::Scalar alpha_in = NAN;
171 : Set::Scalar alpha_out = NAN;
172 : Type type = Type::XYZ;
173 : Unit unit = Unit::Less();
174 : };
175 : }
176 : #endif
|
