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