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