Line data Source code
1 : // If :code:`number_of_inclusions` is specified, then multiple ellipses are specified.
2 : // In this case, each parameter must have number_of_inclusion*M values, where M is the
3 : // number of values specified for the single ellipse case.
4 : //
5 : // .. WARNING::
6 : //
7 : // This class is redundant with :ref:`IC::Ellipsoid` and :ref:`IC::Expression` and will
8 : // probably be consolidated.
9 : //
10 :
11 :
12 : #ifndef IC_ELLIPSE_H_
13 : #define IC_ELLIPSE_H_
14 :
15 : #include "Set/Set.H"
16 : #include "IC/IC.H"
17 : #include "IO/ParmParse.H"
18 :
19 : namespace IC
20 : {
21 : class Ellipse : public IC
22 : {
23 : public:
24 : static constexpr const char* name = "ellipse";
25 :
26 : enum Mollifier {Dirac, Gaussian};
27 :
28 5 : Ellipse (amrex::Vector<amrex::Geometry> &_geom) : IC(_geom) {}
29 5 : Ellipse (amrex::Vector<amrex::Geometry> &_geom, IO::ParmParse &pp, std::string name) : Ellipse(_geom)
30 5 : {pp_queryclass(name,*this);}
31 :
32 52 : void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
33 : {
34 52 : Set::Vector DX(geom[lev].CellSize());
35 52 : amrex::IndexType type = a_field[lev]->ixType();
36 52 : int ncomp = a_field[lev]->nComp();
37 :
38 313 : for (amrex::MFIter mfi(*a_field[lev],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
39 : {
40 261 : amrex::Box bx;
41 522 : if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
42 522 : if (type == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
43 :
44 261 : amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
45 261 : amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
46 :
47 128591 : Set::Vector x;
48 : // NODE
49 257182 : if (type == amrex::IndexType::TheNodeType())
50 : {
51 77231 : AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i)) * geom[lev].CellSize()[0];,
52 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j)) * geom[lev].CellSize()[1];,
53 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k)) * geom[lev].CellSize()[2];);
54 : }
55 102720 : else if (type == amrex::IndexType::TheCellType())
56 : {
57 51360 : AMREX_D_TERM(x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
58 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
59 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
60 : }
61 :
62 128591 : if(number_of_inclusions == 0)
63 : {
64 128591 : Set::Scalar norm = (A[0]*(x-center[0])).lpNorm<2>();
65 128591 : field(i,j,k,0) = 0.5 - 0.5*std::erf(((x-center[0]).transpose() * A[0] * (x-center[0]) - 1.0) / eps[0] / norm);
66 231311 : if (invert) field(i,j,k,0) = 1.0 - field(i,j,k,0);
67 283053 : if (ncomp > 1) field(i,j,k,1) = 1.0 - field(i,j,k,0);
68 : }
69 : else
70 : {
71 : // Base matrix is 0, rest of the inclusions are numbered from 1.
72 0 : Set::Scalar value = 0.0;
73 0 : for (int m = 0; m < number_of_inclusions; m++)
74 : {
75 0 : Set::Scalar norm = (A[m]*(x-center[m])).lpNorm<2>();
76 0 : value += 0.5 - 0.5*std::erf(((x-center[m]).transpose() * A[m] * (x-center[m]) - 1.0) / eps[m] / norm);
77 : // if (field(i,j,k,m+1) < 0.) field(i,j,k,m+1) = 0.;
78 : // if (field(i,j,k,m+1) > 1.) field(i,j,k,m+1) = 1.;
79 : // value += field(i,j,k,m+1);
80 : }
81 0 : field (i,j,k,invert) = 1.0 - value;
82 0 : field (i,j,k,1-invert) = value;
83 0 : if (field(i,j,k,invert) < 0.) field(i,j,k,invert) = 0.;
84 0 : if (field(i,j,k,invert) > 1.) field(i,j,k,invert) = 1.;
85 0 : if (field(i,j,k,1-invert) < 0.) field(i,j,k,1-invert) = 0.;
86 0 : if (field(i,j,k,1-invert) > 1.) field(i,j,k,1-invert) = 1.;
87 : }
88 :
89 128591 : });
90 : }
91 52 : a_field[lev]->FillBoundary();
92 52 : }
93 :
94 : private:
95 : int number_of_inclusions = 1;
96 : amrex::Vector<Set::Vector> center;
97 : amrex::Vector<Set::Matrix> A;
98 : amrex::Vector<Set::Scalar> eps;
99 : int invert = 0;
100 :
101 : public:
102 5 : static void Parse(Ellipse & value, IO::ParmParse & pp)
103 : {
104 10 : amrex::Vector<Set::Scalar> x0;
105 5 : if(!(pp.contains("number_of_inclusions")))
106 : {
107 5 : value.center.resize(0);
108 5 : value.A.resize(0);
109 5 : value.eps.resize(0);
110 :
111 5 : value.number_of_inclusions = 0;
112 5 : pp_queryarr("x0",x0); // Coorinates of ellipse center
113 5 : value.center.push_back(Set::Vector(AMREX_D_DECL(x0[0],x0[1],x0[2])));
114 :
115 : Set::Scalar _eps;
116 5 : pp_query_default("eps",_eps,0.0); // Diffuse boundary thickness
117 5 : value.eps.push_back(_eps);
118 :
119 5 : Set::Matrix _A = Set::Matrix::Zero();
120 5 : if (pp.contains("A"))
121 : {
122 0 : pp_queryarr("A",_A); // DxD square matrix defining an ellipse.
123 0 : value.A.push_back(_A);
124 : }
125 5 : else if (pp.contains("a"))
126 : {
127 5 : Set::Matrix _A = Set::Matrix::Zero();
128 5 : Set::Vector a = Set::Vector::Ones();
129 5 : pp_queryarr("a",a); // If :code:`A` is not defined, then assume a sphere with radius :code:`a`
130 15 : for (int d = 0; d < AMREX_SPACEDIM; d++) _A(d,d) = 1./a(d)/a(d);
131 5 : value.A.push_back(_A);
132 : }
133 : }
134 : else
135 : {
136 0 : pp_query("number_of_inclusions",value.number_of_inclusions); // Number of ellipses
137 0 : if(value.number_of_inclusions < 1) Util::Abort(INFO, "number of inclusions have to be at least 1");
138 :
139 0 : value.center.resize(0);
140 0 : value.A.resize(0);
141 0 : value.eps.resize(0);
142 :
143 0 : if (pp.contains("center") && pp.contains("x0")) Util::Abort(INFO,"Cannot specify both center (depricated) and x0");
144 0 : pp_queryarr("center", x0); // center of the ellipse
145 0 : pp_queryarr("x0", x0); // center of the ellipse
146 :
147 0 : if(x0.size() != value.number_of_inclusions*AMREX_SPACEDIM){
148 0 : Util::Message(INFO, value.number_of_inclusions*AMREX_SPACEDIM);
149 0 : Util::Message(INFO, x0.size());
150 0 : Util::Abort(INFO, "Need centers for all the inclusions");}
151 0 : for (int i = 0; i < x0.size(); i+= AMREX_SPACEDIM)
152 0 : value.center.push_back(Set::Vector(AMREX_D_DECL(x0[i],x0[i+1],x0[i+2])));
153 :
154 0 : if (pp.contains("A"))
155 : {
156 0 : amrex::Vector<Set::Scalar> _A;
157 0 : pp_queryarr("A", _A); // either a vector containing ellipse radii, or a matrix defining the ellipse
158 0 : if(_A.size() != value.number_of_inclusions*AMREX_SPACEDIM*AMREX_SPACEDIM && _A.size() != AMREX_SPACEDIM*AMREX_SPACEDIM)
159 0 : Util::Abort(INFO, "Invalid value of A for ellipse initialization");
160 0 : if(_A.size() == AMREX_SPACEDIM*AMREX_SPACEDIM)
161 : {
162 0 : Set::Matrix _A1 = Set::Matrix::Zero();
163 0 : pp_queryarr("A",_A1); // Same
164 0 : for (int i = 0; i < value.number_of_inclusions; i++)
165 0 : value.A.push_back(_A1);
166 : }
167 : else
168 : {
169 0 : Set::Matrix _A1 = Set::Matrix::Zero();
170 0 : for (int i = 0; i < value.number_of_inclusions; i+= AMREX_SPACEDIM*AMREX_SPACEDIM)
171 : {
172 0 : AMREX_D_PICK( _A1(0,0) = _A[i];
173 : ,
174 : _A1(0,0) = _A[i]; _A1(0,1) = _A[i+1];
175 : _A1(1,0) = _A[i+2]; _A1(1,1) = _A[i+3];
176 : ,
177 : _A1(0,0) = _A[i+0]; _A1(0,1) = _A[i+1]; _A1(0,2) = _A[i+2];
178 : _A1(1,0) = _A[i+3]; _A1(1,1) = _A[i+4]; _A1(1,2) = _A[i+5];
179 : _A1(2,0) = _A[i+6]; _A1(2,1) = _A[i+7]; _A1(2,2) = _A[i+8];
180 : );
181 0 : value.A.push_back(_A1);
182 : }
183 : }
184 : }
185 0 : else if (pp.contains("radius"))
186 : {
187 0 : amrex::Vector<Set::Scalar> _r;
188 0 : pp_queryarr("radius", _r); // Array of radii [depricated]
189 :
190 0 : if(_r.size() != value.number_of_inclusions*AMREX_SPACEDIM && _r.size() != AMREX_SPACEDIM)
191 0 : Util::Abort(INFO, "Invalid value of radius for ellipse initialization");
192 0 : if(_r.size() == AMREX_SPACEDIM)
193 : {
194 0 : Set::Matrix _A1 = Set::Matrix::Zero();
195 0 : for (int i = 0; i< AMREX_SPACEDIM; i++) _A1(i,i) = 1.0/(_r[i]*_r[i]);
196 0 : for (int i = 0; i < value.number_of_inclusions; i++)
197 0 : value.A.push_back(_A1);
198 : }
199 : else
200 : {
201 0 : for (int i = 0; i < value.number_of_inclusions; i++)
202 : {
203 0 : Set::Matrix _A1 = Set::Matrix::Zero();
204 0 : for (int d = 0; d < AMREX_SPACEDIM; d++) _A1(d,d) = 1.0/_r[i+d]/_r[i+d];
205 0 : value.A.push_back(_A1);
206 : }
207 : }
208 : }
209 :
210 0 : amrex::Vector<Set::Scalar> _eps;
211 0 : pp_queryarr("eps",_eps); // Regularization for smooth boundary
212 0 : if(_eps.size() != 1 && _eps.size() != value.number_of_inclusions)
213 0 : Util::Abort(INFO, "Incorrect eps. Check the number of values specified");
214 :
215 0 : if (_eps.size() == 1)
216 0 : for (int i = 0; i < value.number_of_inclusions; i++) value.eps.push_back(_eps[0]);
217 :
218 : else
219 0 : for (int i = 0; i < value.number_of_inclusions; i++) value.eps.push_back(_eps[i]);
220 : }
221 : // Flip the inside and the outside
222 5 : pp_query("invert",value.invert);
223 5 : }
224 : };
225 : }
226 : #endif
|
