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