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