Line data Source code
1 : //
2 : // The Set namespace defines the Alamo datatypes.
3 : // This includes scalars (:code:`Set::Scalar`), vectors (:code:`Set::Vector`), and nD tensors (:code:`Set::Tensor`).
4 : // Scalars are an alias for AMReX :code:`Real`.
5 : // Vectors and tensors are an alias for the Tuxfamily Eigen vectors and matrices, and so you can use
6 : // Eigen operations on them (like :code:`A.determinant()`, etc.)
7 : // See the Eigen documentation for more.
8 : //
9 : // The Set namespace defines the :code:`Matrix4` datatype for fourth order tensors.
10 : // These can have isotropic, major, minor, majorminor, etc., symmetries.
11 : //
12 : // The Set namespace also defines the :code:`Field` datatype, which is where most data is stored.
13 : // It is templated so that :code:`Set::Field<Set::Scalar>` is a field of scalars, :code:`Set::Field<Set::Vector>`
14 : // is a vector field, etc.
15 : // One can also have, for instance, a field of models, like :code:`Set::Field<Model::Solid::Linear::Iostropic>`.
16 : //
17 :
18 : #ifndef SET_SET_
19 : #define SET_SET_
20 :
21 : #include <iomanip>
22 :
23 : #include "Util/Util.H"
24 :
25 : #include "Set/Base.H"
26 : /// \brief A collection of data types and symmetry-reduced data structures
27 :
28 : namespace Set
29 : {
30 :
31 : enum Hypercube {
32 : Node = 0, Edge = 1, Face = 2, Cell = 3
33 : };
34 :
35 : using HC = Hypercube;
36 :
37 : //class Field : public amrex::MultiFab
38 : //{
39 :
40 : //};
41 : template <class T>
42 : class Field : public amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<T>>>>
43 : {
44 : public:
45 1574 : Field() {}
46 : //~Field() {}
47 : Field(int a_levs, const amrex::Vector<amrex::BoxArray> & a_grids, const amrex::Vector<amrex::DistributionMapping> & a_dmap, int a_ncomp, int a_nghost)
48 : {
49 : Define(a_levs,a_grids,a_dmap,a_ncomp,a_nghost);
50 : }
51 : Field(int size) : amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<T>>>>(size) {}
52 : void Define(int a_levs, const amrex::Vector<amrex::BoxArray> & a_grids, const amrex::Vector<amrex::DistributionMapping> & a_dmap, int a_ncomp, int a_nghost)
53 : {
54 : Util::Assert(INFO,TEST(a_levs == a_grids.size()));
55 : Util::Assert(INFO,TEST(a_levs == a_dmap.size()));
56 : this->resize(a_levs);
57 : for (int lev = 0; lev < a_levs; lev++)
58 : {
59 : (*this)[lev].reset(new amrex::FabArray<amrex::BaseFab<T>>(a_grids[lev],a_dmap[lev],a_ncomp,a_nghost));
60 : }
61 : }
62 944 : void Define(int a_lev, const amrex::BoxArray & a_grid, const amrex::DistributionMapping & a_dmap, int a_ncomp, int a_nghost)
63 : {
64 1888 : Util::Assert(INFO,TEST(a_lev < this->size()));
65 944 : (*this)[a_lev].reset(new amrex::FabArray<amrex::BaseFab<T>>(a_grid,a_dmap,a_ncomp,a_nghost));
66 944 : }
67 : int finest_level = 0;
68 199 : void Copy(int /*a_lev*/, amrex::MultiFab &/*a_dst*/, int /*a_dstcomp*/, int /*a_nghost*/) const {}
69 : void CopyFrom(int /*a_lev*/, amrex::MultiFab &/*a_dst*/, int /*a_dstcomp*/, int /*a_nghost*/) const {}
70 : void Add(int /*a_lev*/, amrex::MultiFab &/*a_dst*/, int /*a_dstcomp*/, int /*a_nghost*/) const {}
71 : void AddFrom(int /*a_lev*/, amrex::MultiFab &/*a_src*/, int /*a_srccomp*/, int /*a_nghost*/) const {}
72 311 : int NComp() const {return 0;}
73 0 : std::string Name(int) const {return name;}
74 : std::string name;
75 :
76 0 : amrex::Array4<T> Patch (int lev, amrex::MFIter &mfi) const &
77 : {
78 0 : if (this->size()) return (*(*this)[lev]).array(mfi);
79 0 : else return empty;
80 : }
81 :
82 :
83 : amrex::Array4<T> empty;
84 : };
85 :
86 :
87 : template <typename T>
88 : using Patch = amrex::Array4<T> const&;
89 :
90 :
91 : template<>
92 : ALAMO_SINGLE_DEFINITION
93 1085 : void Field<Set::Vector>::Copy(int a_lev, amrex::MultiFab &a_dst, int a_dstcomp, int a_nghost) const
94 : {
95 3869 : for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
96 : {
97 2784 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
98 2784 : if (bx.ok())
99 : {
100 2784 : amrex::Array4<const Set::Vector> const & src = ((*this)[a_lev])->array(mfi);
101 2784 : amrex::Array4<Set::Scalar> const & dst = a_dst.array(mfi);
102 8718 : for (int n = 0; n < AMREX_SPACEDIM; n++)
103 : {
104 5934 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
105 5936440 : dst(i,j,k,a_dstcomp + n) = src(i,j,k)(n);
106 2968220 : });
107 : }
108 : }
109 : }
110 1085 : }
111 : template<>
112 : ALAMO_SINGLE_DEFINITION
113 14 : void Field<Set::Vector>::CopyFrom(int a_lev, amrex::MultiFab &a_src, int a_dstcomp, int a_nghost) const
114 : {
115 98 : for (amrex::MFIter mfi(a_src, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
116 : {
117 84 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
118 84 : if (bx.ok())
119 : {
120 84 : amrex::Array4<Set::Vector> const & dst = ((*this)[a_lev])->array(mfi);
121 84 : amrex::Array4<const Set::Scalar> const & src = a_src.array(mfi);
122 252 : for (int n = 0; n < AMREX_SPACEDIM; n++)
123 : {
124 168 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
125 287928 : dst(i,j,k)(n) = src(i,j,k,a_dstcomp + n);
126 95976 : });
127 : }
128 : }
129 : }
130 14 : }
131 : template<>
132 : ALAMO_SINGLE_DEFINITION
133 0 : void Field<Set::Vector>::Add(int a_lev, amrex::MultiFab &a_dst, int a_dstcomp, int a_nghost) const
134 : {
135 0 : for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
136 : {
137 0 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
138 0 : if (bx.ok())
139 : {
140 0 : amrex::Array4<const Set::Vector> const & src = ((*this)[a_lev])->array(mfi);
141 0 : amrex::Array4<Set::Scalar> const & dst = a_dst.array(mfi);
142 0 : for (int n = 0; n < AMREX_SPACEDIM; n++)
143 : {
144 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
145 0 : dst(i,j,k,a_dstcomp + n) += src(i,j,k)(n);
146 0 : });
147 : }
148 : }
149 : }
150 0 : }
151 : template<>
152 : ALAMO_SINGLE_DEFINITION
153 1372 : void Field<Set::Vector>::AddFrom(int a_lev, amrex::MultiFab &a_src, int a_srccomp, int a_nghost) const
154 : {
155 2877 : for (amrex::MFIter mfi(a_src, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
156 : {
157 1505 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
158 1505 : if (bx.ok())
159 : {
160 1505 : amrex::Array4<Set::Vector> const & dst = ((*this)[a_lev])->array(mfi);
161 1505 : amrex::Array4<const Set::Scalar> const & src = a_src.array(mfi);
162 4815 : for (int n = 0; n < AMREX_SPACEDIM; n++)
163 : {
164 3310 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
165 3373070 : dst(i,j,k)(n) += src(i,j,k,a_srccomp + n);
166 1124358 : });
167 : }
168 : }
169 : }
170 1372 : }
171 :
172 2295 : template<> ALAMO_SINGLE_DEFINITION int Field<Set::Vector>::NComp() const {return AMREX_SPACEDIM;}
173 416 : template<> ALAMO_SINGLE_DEFINITION std::string Field<Set::Vector>::Name(int i) const
174 : {
175 : #if AMREX_SPACEDIM>0
176 416 : if (i==0) return name + "_x";
177 : #endif
178 : #if AMREX_SPACEDIM>1
179 241 : else if (i==1) return name + "_y";
180 : #endif
181 : #if AMREX_SPACEDIM>2
182 66 : else if (i==2) return name + "_z";
183 : #endif
184 0 : else Util::Abort(INFO,"Invalid component");
185 0 : return "";
186 : }
187 :
188 : template<>
189 : ALAMO_SINGLE_DEFINITION
190 532 : void Field<Set::Matrix>::Copy(int a_lev, amrex::MultiFab &a_dst, int a_dstcomp, int a_nghost) const
191 : {
192 2186 : for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
193 : {
194 1654 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
195 1654 : if (bx.ok())
196 : {
197 1654 : amrex::Array4<const Set::Matrix> const & src = ((*this)[a_lev])->array(mfi);
198 1654 : amrex::Array4<Set::Scalar> const & dst = a_dst.array(mfi);
199 8600 : for (int n = 0; n < AMREX_SPACEDIM*AMREX_SPACEDIM; n++)
200 : {
201 6946 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
202 6224180 : dst(i,j,k,a_dstcomp + n) = src(i,j,k)(n/AMREX_SPACEDIM,n%AMREX_SPACEDIM);
203 3112090 : });
204 : }
205 : }
206 : }
207 532 : }
208 : template<>
209 1870 : ALAMO_SINGLE_DEFINITION int Field<Set::Matrix>::NComp() const {return AMREX_SPACEDIM*AMREX_SPACEDIM;}
210 1002 : template<> ALAMO_SINGLE_DEFINITION std::string Field<Set::Matrix>::Name(int i) const
211 : {
212 : #if AMREX_SPACEDIM==1
213 : if (i==0) return name + "_xx";
214 : #elif AMREX_SPACEDIM==2
215 408 : if (i==0) return name + "_xx";
216 306 : else if (i==1) return name + "_xy";
217 204 : else if (i==2) return name + "_yx";
218 102 : else if (i==3) return name + "_yy";
219 : #elif AMREX_SPACEDIM==3
220 594 : if (i==0) return name + "_xx";
221 528 : else if (i==1) return name + "_xy";
222 462 : else if (i==2) return name + "_xz";
223 396 : else if (i==3) return name + "_yx";
224 330 : else if (i==4) return name + "_yy";
225 264 : else if (i==5) return name + "_yz";
226 198 : else if (i==6) return name + "_zx";
227 132 : else if (i==7) return name + "_zy";
228 66 : else if (i==8) return name + "_zz";
229 : #endif
230 0 : else Util::Abort(INFO, "Invalid component");
231 0 : return "";
232 : }
233 :
234 :
235 : template <>
236 : class Field<Set::Scalar> : public amrex::Vector<std::unique_ptr<amrex::MultiFab>>
237 : {
238 : public:
239 1460 : Field() {}
240 12 : Field(int size) : amrex::Vector<std::unique_ptr<amrex::MultiFab>>(size) {}
241 : void Define(int a_levs, const amrex::Vector<amrex::BoxArray> & a_grids, const amrex::Vector<amrex::DistributionMapping> & a_dmap, int a_ncomp, int a_nghost)
242 : {
243 : Util::Assert(INFO,TEST(a_levs == a_grids.size()));
244 : Util::Assert(INFO,TEST(a_levs == a_dmap.size()));
245 : this->resize(a_levs);
246 : for (int lev = 0; lev < a_levs; lev++)
247 : {
248 : (*this)[lev].reset(new amrex::MultiFab(a_grids[lev],a_dmap[lev],a_ncomp,a_nghost));
249 : }
250 : }
251 992 : void Define(int a_lev, const amrex::BoxArray & a_grid, const amrex::DistributionMapping & a_dmap, int a_ncomp, int a_nghost)
252 : {
253 1984 : Util::Assert(INFO,TEST(a_lev < this->size()));
254 992 : (*this)[a_lev].reset(new amrex::MultiFab(a_grid,a_dmap,a_ncomp,a_nghost));
255 992 : }
256 : int finest_level = 0;
257 :
258 : void Copy(int /*a_lev*/, amrex::MultiFab &/*a_dst*/, int /*a_dstcomp*/, int /*a_nghost*/) const
259 : {Util::Abort(INFO,"This should never get called");}
260 : int NComp() const
261 : {Util::Abort(INFO,"This should never be called");return -1;}
262 :
263 323162 : amrex::Array4<Set::Scalar> Patch (int lev, amrex::MFIter &mfi) const &
264 : {
265 323162 : if (this->size()) return (*(*this)[lev]).array(mfi);
266 30662 : else return empty;
267 : }
268 :
269 : amrex::Array4<Set::Scalar> empty;
270 : };
271 :
272 :
273 : }
274 :
275 : namespace Util
276 : {
277 : Set::Scalar Random();
278 : Set::Scalar Gaussian(amrex::Real mean,amrex::Real std_deviation);
279 : }
280 :
281 : namespace Set
282 : {
283 : namespace Constant
284 : {
285 : // Pi
286 : static const Set::Scalar Pi = 3.14159265359;
287 : // Gas constant
288 : static const Set::Scalar Rg = 8.31446261815;
289 : }
290 : }
291 :
292 : #include "Set/Matrix3.H"
293 : #include "Set/Matrix4.H"
294 :
295 :
296 : #endif
|
