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