Line data Source code
1 : #ifndef INTEGRATOR_BASEFIELD_H
2 : #define INTEGRATOR_BASEFIELD_H
3 :
4 : #include "AMReX_FillPatchUtil.H"
5 :
6 : #include "Util/Util.H"
7 : #include "Set/Set.H"
8 : #include "Numeric/Interpolator/NodeBilinear.H"
9 :
10 : namespace Integrator
11 : {
12 :
13 : class BaseField
14 : {
15 : public:
16 98 : virtual ~BaseField() = default;
17 :
18 : virtual void RemakeLevel (int lev,
19 : amrex::Real time,
20 : const amrex::BoxArray& cgrids,
21 : const amrex::DistributionMapping& dm) = 0;
22 : virtual void MakeNewLevelFromCoarse (int lev,
23 : amrex::Real time,
24 : const amrex::BoxArray& cgrids,
25 : const amrex::DistributionMapping& dm) = 0;
26 : virtual void MakeNewLevelFromScratch (int lev,
27 : amrex::Real t,
28 : const amrex::BoxArray& cgrids,
29 : const amrex::DistributionMapping& dm) = 0;
30 : virtual void SetFinestLevel(const int a_finestlevel) = 0;
31 : virtual void FillPatch(const int lev, const Set::Scalar time)=0;
32 : virtual void FillBoundary(const int lev, Set::Scalar time)=0;
33 : virtual void AverageDown(const int lev, amrex::IntVect refRatio) = 0;
34 :
35 : virtual int NComp() = 0;
36 : virtual void Copy(int /*a_lev*/, amrex::MultiFab &/*a_dst*/, int /*a_dstcomp*/, int /*a_nghost*/) = 0;
37 : bool writeout = false;
38 : virtual std::string Name(int) = 0;
39 : virtual void setName(std::string a_name) = 0;
40 : bool evolving = true;
41 : virtual void setBC(void * a_bc) = 0;
42 : virtual void * getBC() = 0;
43 : Set::Hypercube m_gridtype;
44 : };
45 :
46 : template<class T>
47 : class Field : public BaseField
48 : {
49 : public:
50 : class EmptyBC
51 : {
52 : public:
53 45349 : virtual void operator () (amrex::FabArray<amrex::BaseFab<T>> & /*mf*/,
54 : int /*dcomp*/, int /*ncomp*/, amrex::IntVect const& /*nghost*/,
55 : amrex::Real /*time*/, int /*bccomp*/)
56 45349 : { /* Do nothing - this is a shell to satisfy the FillPatch functions.*/ }
57 13028 : amrex::BCRec GetBCRec() {return amrex::BCRec();}
58 : };
59 :
60 272 : Field(Set::Field<T> & a_field,
61 : const amrex::Vector<amrex::Geometry> &a_geom,
62 : const amrex::Vector<amrex::IntVect> &a_refRatio,
63 : int a_ncomp, int a_nghost) :
64 272 : m_field(a_field), m_geom(a_geom), m_refRatio(a_refRatio),
65 272 : m_ncomp(a_ncomp), m_nghost(a_nghost)
66 272 : {}
67 :
68 : void
69 42051 : FillPatch (int lev, amrex::Real time,
70 : amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<T>>>> &source_mf,
71 : amrex::FabArray<amrex::BaseFab<T>> &destination_mf,
72 : int icomp)
73 : {
74 42051 : EmptyBC physbc;
75 42051 : if (lev == 0)
76 : {
77 38836 : amrex::Vector<amrex::FabArray<amrex::BaseFab<T>>*> smf;
78 38836 : smf.push_back(source_mf[lev].get());
79 38836 : amrex::Vector<amrex::Real> stime;
80 38836 : stime.push_back(time);
81 :
82 38836 : if (m_gridtype == Set::Hypercube::Cell)
83 : {
84 0 : m_bc->define(m_geom[lev]);
85 0 : amrex::FillPatchSingleLevel(destination_mf, time, smf, stime,
86 0 : 0, icomp, destination_mf.nComp(), m_geom[lev],
87 0 : *m_bc, 0);
88 : }
89 38836 : else if (m_gridtype == Set::Hypercube::Node)
90 : {
91 38836 : amrex::FillPatchSingleLevel(destination_mf, time, smf, stime,
92 38836 : 0, icomp, destination_mf.nComp(), m_geom[lev],
93 : physbc, 0);
94 : }
95 38836 : }
96 : else
97 : {
98 3215 : amrex::Vector<amrex::FabArray<amrex::BaseFab<T>>*> cmf, fmf;
99 3215 : cmf.push_back(source_mf[lev-1].get());
100 3215 : fmf.push_back(source_mf[lev].get());
101 3215 : amrex::Vector<amrex::Real> ctime, ftime;
102 3215 : ctime.push_back(time);
103 3215 : ftime.push_back(time);
104 :
105 3215 : amrex::Vector<amrex::BCRec> bcs(destination_mf.nComp(), physbc.GetBCRec());
106 :
107 6430 : if (destination_mf.boxArray().ixType() == amrex::IndexType::TheNodeType())
108 : {
109 3215 : amrex::FillPatchTwoLevels(destination_mf, time, cmf, ctime, fmf, ftime,
110 3215 : 0, icomp, destination_mf.nComp(), m_geom[lev-1], m_geom[lev],
111 : physbc, 0,
112 : physbc, 0,
113 3215 : m_refRatio[lev-1],
114 : &node_bilinear, bcs, 0);
115 :
116 : }
117 : else
118 : {
119 0 : Util::Abort(INFO,"This implementation is not yet complete.");
120 : /*
121 : Numeric::Interpolator::CellBilinear<T> *mapper
122 : = new Numeric::Interpolator::CellBilinear<T>();
123 : m_bc->define(m_geom[lev]);
124 : amrex::FillPatchTwoLevels(destination_mf, time, cmf, ctime, fmf, ftime,
125 : 0, icomp, destination_mf.nComp(), m_geom[lev-1], m_geom[lev],
126 : *m_bc, 0,
127 : *m_bc, 0,
128 : m_refRatio[lev-1],
129 : mapper, bcs, 0);
130 : */
131 : }
132 3214 : }
133 :
134 42050 : }
135 : void
136 41533 : FillPatch (int lev, amrex::Real time) override
137 : {
138 41533 : this->FillPatch(lev,time,m_field,*(m_field[lev]),0);
139 41532 : }
140 : void
141 42 : FillCoarsePatch (int lev,
142 : amrex::Real time,
143 : int icomp,
144 : int ncomp)
145 : {
146 : BL_PROFILE("Integrator::FillCoarsePatch");
147 : AMREX_ASSERT(lev > 0);
148 :
149 42 : if (m_gridtype == Set::Hypercube::Node)
150 : {
151 336 : Util::Assert(INFO, TEST(m_field[lev]->boxArray().ixType() == amrex::IndexType::TheNodeType()));
152 42 : EmptyBC physbc;
153 42 : amrex::Vector<amrex::FabArray<amrex::BaseFab<T>> *> cmf;
154 42 : cmf.push_back(m_field[lev-1].get());
155 42 : amrex::Vector<amrex::Real> ctime;
156 42 : ctime.push_back(time);
157 42 : amrex::Vector<amrex::BCRec> bcs(ncomp, physbc.GetBCRec());
158 42 : amrex::InterpFromCoarseLevel(*m_field[lev], time, *cmf[0], 0, icomp, ncomp,
159 42 : m_geom[lev-1], m_geom[lev],
160 : physbc, 0,
161 : physbc, 0,
162 42 : m_refRatio[lev-1],
163 : &node_bilinear, bcs, 0);
164 42 : }
165 0 : else if (m_gridtype == Set::Hypercube::Cell)
166 : {
167 0 : Util::Abort(INFO,"This implementation is not yet complete");
168 : /*
169 : amrex::Vector<amrex::FabArray<amrex::BaseFab<T>> *> cmf;
170 : cmf.push_back(m_field[lev-1].get());
171 : amrex::Vector<amrex::Real> ctime;
172 : ctime.push_back(time);
173 : m_bc->define(m_geom[lev]);
174 : Numeric::Interpolator::CellBilinear<T> *mapper = new Numeric::Interpolator::CellBilinear<T>();
175 : amrex::Vector<amrex::BCRec> bcs(ncomp, m_bc->GetBCRec());
176 : amrex::InterpFromCoarseLevel(*m_field[lev], time, *cmf[0], 0, icomp, ncomp,
177 : m_geom[lev-1], m_geom[lev],
178 : *m_bc, 0,
179 : *m_bc, 0,
180 : m_refRatio[lev-1],
181 : mapper, bcs, 0);
182 : */
183 : }
184 : else
185 0 : Util::Abort(INFO,"Invalid grid");
186 42 : }
187 :
188 :
189 518 : virtual void RemakeLevel (int lev,
190 : amrex::Real time,
191 : const amrex::BoxArray& cgrids,
192 : const amrex::DistributionMapping& dm) override
193 : {
194 518 : if (m_gridtype == Set::Hypercube::Node)
195 : {
196 4144 : Util::Assert(INFO, TEST(m_field[lev]->boxArray().ixType() == amrex::IndexType::TheNodeType()));
197 518 : amrex::BoxArray ngrids = cgrids;
198 518 : ngrids.convert(amrex::IntVect::TheNodeVector());
199 518 : amrex::FabArray<amrex::BaseFab<T>> new_state(ngrids, dm, m_ncomp, m_nghost);
200 518 : this->FillPatch(lev, time, m_field, new_state, 0);
201 518 : std::swap(new_state, *m_field[lev]);
202 518 : }
203 0 : else if (m_gridtype == Set::Hypercube::Cell)
204 : {
205 0 : amrex::FabArray<amrex::BaseFab<T>> new_state(cgrids, dm, m_ncomp, m_nghost);
206 0 : this->FillPatch(lev, time, m_field, new_state, 0);
207 0 : std::swap(new_state, *m_field[lev]);
208 0 : }
209 : else
210 0 : Util::Abort(INFO,"Invalid grid");
211 518 : }
212 :
213 42 : virtual void MakeNewLevelFromCoarse (int lev,
214 : amrex::Real time,
215 : const amrex::BoxArray& cgrids,
216 : const amrex::DistributionMapping& dm) override
217 : {
218 42 : if (m_gridtype == Set::Hypercube::Node)
219 : {
220 42 : amrex::BoxArray ngrids = cgrids;
221 42 : ngrids.convert(amrex::IntVect::TheNodeVector());
222 42 : m_field[lev].reset(new amrex::FabArray<amrex::BaseFab<T>>(ngrids,dm,m_ncomp,m_nghost));
223 42 : FillCoarsePatch(lev,time,0,m_ncomp);
224 42 : }
225 0 : else if (m_gridtype == Set::Hypercube::Cell)
226 : {
227 0 : m_field[lev].reset(new amrex::FabArray<amrex::BaseFab<T>>(cgrids,dm,m_ncomp,m_nghost));
228 0 : FillCoarsePatch(lev,time,0,m_ncomp);
229 : }
230 : else
231 0 : Util::Abort(INFO,"Invalid grid");
232 42 : }
233 :
234 1417 : virtual void MakeNewLevelFromScratch (int lev,
235 : amrex::Real /*time*/,
236 : const amrex::BoxArray& cgrids,
237 : const amrex::DistributionMapping& dm) override
238 : {
239 1417 : if (m_gridtype == Set::Hypercube::Node)
240 : {
241 1417 : amrex::BoxArray ngrids = cgrids;
242 1417 : ngrids.convert(amrex::IntVect::TheNodeVector());
243 1417 : m_field[lev].reset(new amrex::FabArray<amrex::BaseFab<T>>(ngrids,dm,m_ncomp,m_nghost));
244 1417 : m_field[lev]->setVal(T::Zero());
245 1417 : }
246 0 : else if (m_gridtype == Set::Hypercube::Cell)
247 : {
248 0 : m_field[lev].reset(new amrex::FabArray<amrex::BaseFab<T>>(cgrids,dm,m_ncomp,m_nghost));
249 0 : m_field[lev]->setVal(T::Zero());
250 : }
251 : else
252 0 : Util::Abort(INFO,"Invalid grid");
253 1417 : }
254 :
255 125581 : virtual void SetFinestLevel(const int a_finestlevel) override
256 : {
257 125581 : m_field.finest_level = a_finestlevel;
258 125581 : }
259 :
260 1417 : virtual void FillBoundary(const int lev, Set::Scalar time) override
261 : {
262 1417 : Util::RealFillBoundary(*m_field[lev],m_geom[lev]);
263 1417 : if (m_gridtype == Set::Hypercube::Cell)
264 : {
265 0 : m_bc->define(m_geom[lev]);
266 0 : for (amrex::MFIter mfi(*m_field[lev],true); mfi.isValid(); ++mfi)
267 : {
268 0 : const amrex::Box& box = mfi.tilebox();
269 0 : amrex::BaseFab<T> &patch = (*(m_field[lev]))[mfi];
270 0 : m_bc->FillBoundary(patch,box,m_nghost,0,m_ncomp,time);
271 : }
272 : //
273 : }
274 1417 : }
275 :
276 1336 : virtual void AverageDown(const int lev, amrex::IntVect refRatio) override
277 : {
278 1336 : if (m_gridtype == Set::Hypercube::Node)
279 1336 : amrex::average_down_nodal(*m_field[lev+1],*m_field[lev],refRatio);
280 0 : else if (m_gridtype == Set::Hypercube::Cell)
281 0 : Util::Abort(INFO,"This is not yet working");
282 : //amrex::average_down(*m_field[lev+1],*m_field[lev],0,m_ncomp,refRatio);
283 : //amrex::average_down(*m_field[lev+1],*m_field[lev],0,m_ncomp,refRatio);
284 : else
285 0 : Util::Abort(INFO,"Incorrect grid");
286 1336 : }
287 :
288 : private:
289 :
290 : Set::Field<T> &m_field;
291 : const amrex::Vector<amrex::Geometry> &m_geom;
292 : const amrex::Vector<amrex::IntVect> &m_refRatio;
293 : const int m_ncomp, m_nghost;
294 : BC::BC<T> *m_bc;
295 :
296 : Numeric::Interpolator::NodeBilinear<T> node_bilinear;
297 :
298 : public:
299 10055 : virtual int NComp() override {
300 10055 : return m_field.NComp();
301 : }
302 3933 : virtual void Copy(int a_lev, amrex::MultiFab &a_dst, int a_dstcomp, int a_nghost) override
303 : {
304 3933 : m_field.Copy(a_lev,a_dst,a_dstcomp,a_nghost);
305 3933 : }
306 272 : virtual void setName(std::string a_name) override {
307 272 : m_field.name = a_name;
308 272 : }
309 3822 : virtual std::string Name(int i) override {
310 3822 : return m_field.Name(i);
311 : }
312 0 : virtual void setBC(void * a_bc) override {
313 0 : m_bc = (BC::BC<T> *)(a_bc);
314 0 : }
315 0 : virtual void * getBC() override {
316 0 : return (void*)m_bc;
317 : }
318 :
319 : };
320 :
321 : }
322 :
323 : #endif
|
