Chombo + EB + MF  3.2
IndexedMomentsImplem.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _INDEXEDMOMENTSIMPLEM_H_
12 #define _INDEXEDMOMENTSIMPLEM_H_
13 
14 #include "MomentIterator.H"
15 #include "EB_TYPEDEFS.H"
16 #include "IndexTM.H"
17 #include "Factorial.H"
18 #include "CH_Timer.H"
19 #include "parstream.H"
20 #include <map>
21 #include <utility>
22 #include "NamespaceHeader.H"
23 
24 template <int Dim, int P> bool IndexedMoments<Dim, P>::s_staticsSet = false;
26 template <int Dim, int P> int IndexedMoments<Dim, P>::s_size = 0;
27 using std::endl;
28 
29 ///
30 /**
31  Moments are centered at the center of the cell.
32  For each of these moments I shift them to the lowest
33  corner of the cell, where I know what the bounds of the integrals
34  is (lower bound always zero, upper bound = dx^d dx^px dx^py dx^pz
35  If the shifted moment is out of bounds, I bound it throw a message.
36  tolerance is about verbosity.
37 **/
38 template<int Dim, int ORDER>
39 void
41  Vector<IndexTM<int, Dim> > & a_bogusPowers,
42  const Real & a_dx,
43  const Real & a_tolerance,
44  const bool & a_ebMoment,
45  const bool & a_bindMoments)
46 {
48  //eb moments can mysteriously get the wrong signs
49  if(a_ebMoment && a_bindMoments)
50  {
51  if(a_moments[IndexTM<int,Dim>::Zero] < 0.)
52  {
53  a_moments *= -1.0;
54  }
55  }
56  a_bogusPowers.resize(0);
57  Real cellvol = 1;
58  for(int idir = 0; idir < Dim; idir++)
59  {
60  cellvol *= a_dx;
61  }
62 
63  //shift moments to lowest corner
64  IndexTM<Real, Dim> shiftvec;
65  shiftvec.setAll(0.5*a_dx);
66  a_moments.shift(shiftvec);
67 
68  //because they are all shifted, to the lowest
69  //min is always zero.
70  const Real zerval = 0;
71  for(momit.reset(); momit.ok(); ++momit)
72  {
73  const IndexTM<int, Dim>& p = momit();
74  //this makes max = dx^p
75  Real maxval = POW<Dim>(a_dx, p);
76  //because this is a vol integral...
77  maxval *= cellvol;
78  //this is a reference so it will change the input
79  Real& momval = a_moments[p];
80  if(momval < zerval)
81  {
82  if((zerval - momval) > a_tolerance*maxval)
83  {
84  a_bogusPowers.push_back(p);
85  // MayDay::Error("minval break low");
86  }
87  if(a_bindMoments)
88  {
89  momval = zerval;
90  }
91  }
92  else if((momval > maxval) && (!a_ebMoment))
93  {
94  if((momval - maxval) > a_tolerance*maxval)
95  {
96  a_bogusPowers.push_back(p);
97  // MayDay::Error("maxval break");
98  }
99  if(a_bindMoments)
100  {
101  momval = maxval;
102  }
103  }
104  }
105 
106  //shift moments back to highest corner
107  shiftvec.setAll(-a_dx);
108  a_moments.shift(shiftvec);
109 
110  //because the moments are shifted to the highest corner
111  //even powers will always be positive, odd ones will always be negtive
112  for(momit.reset(); momit.ok(); ++momit)
113  {
114  const IndexTM<int, Dim>& p = momit();
115  int psum = p.sum();
116  bool pEven = (psum % 2 == 0);
117  Real& momval = a_moments[p];
118  //this makes max = dx^p
119  Real maxval = POW<Dim>(a_dx, p);
120  //because this is a vol integral...
121  maxval *= cellvol;
122  if(pEven && (momval < zerval) && (!a_ebMoment))
123  {
124  //p has an sum and therefore the moment should be positive
125  if(momval < -a_tolerance*maxval )
126  {
127  a_bogusPowers.push_back(p);
128  // MayDay::Error("high corner break even");
129  }
130  if(a_bindMoments)
131  {
132  momval = 0;
133  }
134  }
135  else if(!pEven && (momval > zerval) && (!a_ebMoment))
136  {
137  //p has an odd sum and therefore the moment should be negative
138  if(momval > a_tolerance*maxval)
139  {
140  a_bogusPowers.push_back(p);
141  // MayDay::Error("high corner break odd");
142  }
143  if(a_bindMoments)
144  {
145  momval = 0;
146  }
147  }
148  }
149 
150  //shift moments back to cell center
151  shiftvec.setAll(0.5*a_dx);
152  a_moments.shift(shiftvec);
153 }
154 
155 template <int Dim, int P>
157  {
158  {1, 2, 3, 4, 5, 6, 7, 8, 9}, // D=1
159  {1, 3, 6, 10, 15, 21, 28, 36, 45}, // D=2
160  {1, 4, 10, 20, 35, 56, 84, 120, 165} // D=3
161  };
162 using std::endl;
163 /////
164 template <const int Dim, const int P>
165 Real
168  const map<IndexTM<int, Dim>, Real>& a_mapin) const
169 {
170  /**/
171  typename std::map<IndexTM<int, Dim>, Real>::const_iterator iterinator;
172  iterinator = a_mapin.find(a_mono);
173 
174  Real moment = 1.23456789e10;
175 
176  if (iterinator != a_mapin.end())
177  {
178  moment = iterinator->second;
179 
180  }
181  else
182  {
183  pout() << "monomial = " ;
184  for(int idir = 0; idir < Dim; idir++)
185  {
186  pout() << a_mono[idir] << " ";
187  }
188  pout() << endl;
189  MayDay::Abort("getMoments in IndexedMoments: moment not found in map");
190  }
191  return moment;
192 }
193 
194 ///multiply each component by factor
195 template <int Dim, int P>
198 operator*=(const Real& a_factor)
199 {
200  for(MomentIterator<Dim,P> it; it.ok(); ++it)
201  {
202  (*this)[it()] *= a_factor;
203  }
204  return *this;
205 }
206 ///add compoenentwise
207 template <int Dim, int P>
211 {
212  for(MomentIterator<Dim,P> it; it.ok(); ++it)
213  {
214  (*this)[it()] += a_input[it()];
215  }
216  return *this;
217 }
218 ///shift moment by the input distance.
219 template <int Dim, int P>
220 void
222 shift(const IndexTM<Real, Dim>& a_distance)
223 {
224  CH_TIME("IndexedMoment::shift");
225  IndexTM<Real, Dim> x = a_distance;
226 
227  IndexedMoments<Dim, P> retval = *this;
228  for(MomentIterator<Dim,P> itouter; itouter.ok(); ++itouter)
229  {
230  IndexTM<int, Dim> r = itouter();
231  // - use the binomial shift rule to move them to a_x0
232  Real moment = 0.0;
233  for(MomentIterator<Dim,P> itinner; itinner.ok(); ++itinner)
234  {
235  IndexTM<int, Dim> q = itinner();
236 
237  if(q.componentwiseLE(r))
238  {
239  IndexTM<int, Dim> p = r - q;
240  Real m0 = (*this)[q];
241  Real xpow = 1;
242  for(int idir = 0; idir < Dim; idir++)
243  {
244  xpow *= POW(x[idir], p[idir]);
245  }
246 
247  moment += pCk<Dim>(r, q) * xpow * m0;
248 
249  } //end if (q <= r)
250  }//end loop over inner moments
251  retval[r] = moment;
252  }
253  *this = retval;
254 }
255 
256 /*******/
257 /// for use with irregnode
258 template <int Dim, int P>
261 operator=(const map<IndexTM<int, Dim>, Real>& a_mapin)
262 {
263 
264  for(MomentIterator<Dim, P> iter; iter.ok(); ++iter)
265  {
266  IndexTM<int, Dim> mono = iter();
267  Real moment= getMoment(iter(), a_mapin);
268  (*this)[iter()] = moment;
269  }
270  return *this;
271 }
272 
273 template <int Dim, int P>
274 void
276 spout() const
277 {
278  pout() << "index \t moment " << endl;
279  for(MomentIterator<Dim, P> iter; iter.ok(); ++iter)
280  {
281  pout() << iter() << "\t" << (*this)[iter()] << endl;
282  }
283 }
284 
285 template <int Dim>
286 bool allPositive(const IndexTM<int, Dim>& a_index)
287 {
288  bool allpos = true;
289  for(int idir = 0; idir < Dim; idir++)
290  {
291  if(a_index[idir] < 0)
292  {
293  allpos = false;
294  }
295  }
296  return allpos;
297 }
298 
299 
300 template <int Dim, int P>
301 int
304 {
305  CH_assert(a_index.sum() <= P);
306  CH_assert(allPositive(a_index));
307 
308  // Formula for the layout index with 0th dim first, then 1st, etc.
309  int index = a_index[0];
310  int rem = P;
311  for (int d=Dim-1; d > 0; --d)
312  {
313  index += s_max_sizes[d][rem] - s_max_sizes[d][rem-a_index[d]];
314  rem -= a_index[d];
315  }
316 
317  return index;
318 }
319 /*******/
320 template <int Dim, int P>
321 void
324 {
325  setSize();
326  setMultiIndicies();
327  s_staticsSet= true;
328 }
329 
330 /*******/
331 template <int Dim, int P>
334 {
335  if(!s_staticsSet)
336  {
337  setStatics();
338  }
339 
340  m_moms.resize(s_size);
341  m_isRegular = 0;
342 }
343 
344 /*******/
345 template <int Dim, int P>
346 void
348 setRegular(const Real a_dx)
349 {
350  CH_assert(a_dx > 0);
351  m_isRegular = 1;
352 
353  // Get the set of multi-indexes in the correct order
355 
356  // pout() << "Regular moments:" << endl;
357  for (int ix=0; ix < mix.size(); ++ix)
358  {
359  IndexTM<int,Dim> index = mix[ix];
360 
361  bool even = true;
362  for (int d=0; (even) && (d < Dim); ++d)
363  even = even && !(index[d] % 2);
364 
365  //! Only moments with an even multi-index are non-zero
366  Real moment = 0;
367  if (even)
368  {
369  moment = POW(0.5 * a_dx, index.sum()) * POW(a_dx, Dim);
370  for (int d=0; d < Dim; ++d)
371  moment /= (Real) (index[d]+1);
372  }
373  m_moms[ix] = moment;
374 
375  // pout() << ix << ": " << index << " = " << moment << endl;
376  }
377 
378 }
379 /*******/
380 template <int Dim, int P>
381 void
384 {
385  CH_assert(Dim <= SpaceDim);
386  CH_assert(Dim > 0);
388  CH_assert(P >= 0);
389 
390 
391  s_size = s_max_sizes[Dim-1][P];
392 }
393 /*******/
394 template <int Dim, int P>
395 void
398 {
399  s_multiIndicies.resize(s_size);
400 
402  for (int ix=0; ix < s_size; ++ix)
403  {
404  // If the sum is too large, shift extras to the right
405  for (int d=0; (index.sum() > P) && (d < Dim-1); ++d)
406  {
407  index[d] = 0;
408  ++index[d+1];
409  }
410  s_multiIndicies[ix] = index;
411  ++index[0];
412  }
413 }
414 /****/
415 template <int Dim, int P>
416 void
419  const IndexedMoments<Dim, P> & a_CB)
420 {
421  this->setToZero();
422  //here we are computing sum(m_p) = integral(ni(x-xbar)^p1/p! (x-xbar)^p2)
423  // = sum_p2 ((n^p2)/p2!) m_p+p12
424  for(MomentIterator<Dim, P> momitOuter; momitOuter.ok(); ++momitOuter)
425  {
426  //this is where the answer goes
427  IndexTM<int,SpaceDim> po = momitOuter();
428 
429  for(MomentIterator<Dim, P> momitInner; momitInner.ok(); ++momitInner)
430  {
431  IndexTM<int,SpaceDim> pi = momitInner();
432  IndexTM<int,SpaceDim> psum = po + pi;
433  if(psum.sum() <= P)
434  {
435  Real incr = a_CA[pi]*a_CB[po];
436  (*this)[psum] += incr;
437  }
438  }
439  }
440 
441 }
442 /****/
443 template <int Dim, int P>
444 void
447 {
448  for(MomentIterator<Dim, P> momit; momit.ok(); ++momit)
449  {
450  IndexTM<int,SpaceDim> p = momit();
452  (*this)[p] /= factorial;
453  }
454 }
455 /****/
456 template <int Dim, int P>
457 void
460 {
461  for(MomentIterator<Dim, P> momit; momit.ok(); ++momit)
462  {
463  IndexTM<int,SpaceDim> p = momit();
465  (*this)[p] *= factorial;
466  }
467 }
468 /*******/
469 #include "NamespaceFooter.H"
470 #endif
std::ostream & pout()
Use this in place of std::cout for program output.
bool componentwiseLE(const IndexTM &a_s)
Definition: IndexTM.H:257
static int indexOf(IndexTM< int, Dim > a_index)
Calculate what linear index this multi-index is.
Definition: IndexedMomentsImplem.H:303
Definition: IndexedMoments.H:31
void setRegular(const Real a_dx)
set to a regular IndexTM
Definition: IndexedMomentsImplem.H:348
#define CH_assert(cond)
Definition: CHArray.H:37
void setAll(T a_val)
Definition: IndexTMI.H:194
void shift(const IndexTM< Real, Dim > &a_distance)
shift moment by the input distance.
Definition: IndexedMomentsImplem.H:222
bool allPositive(const IndexTM< int, Dim > &a_index)
Definition: IndexedMomentsImplem.H:286
IndexedMoments()
constructor—make statics first time called
Definition: IndexedMomentsImplem.H:333
one dimensional dynamic array
Definition: Vector.H:53
static const Vector< IndexTM< int, Dim > > & getMonomialPowers()
monomial powers
Definition: IndexedMoments.H:134
Real pfactorial(const IndexTM< int, Dim > p)
Calculates factorials for a multinomial.
Definition: Factorial.H:61
virtual void reset()
Definition: MomentIterator.H:74
Real POW(const Real &a_x, const int &a_p)
computes x^p
Definition: Factorial.H:33
const int SpaceDim
Definition: SPACE.H:38
Definition: IndexTM.H:36
IndexedMoments< Dim, P > & operator*=(const Real &a_factor)
multiply each component by constant
Definition: IndexedMomentsImplem.H:198
void spout() const
outputs the current state to pout() (a la parstream.H)
Definition: IndexedMomentsImplem.H:276
IndexedMoments< Dim, P > & operator=(const map< IndexTM< int, Dim >, Real > &a_mapin)
for use with irregnode
Definition: IndexedMomentsImplem.H:261
static void setStatics()
Definition: IndexedMomentsImplem.H:323
#define CH_TIME(name)
Definition: CH_Timer.H:82
void multiplyByFactorial()
multiply each entry by p!
Definition: IndexedMomentsImplem.H:459
double Real
Definition: REAL.H:33
iterates through the indices of a IndexedMoment
Definition: MomentIterator.H:56
void setToTruncatedMultiply(const IndexedMoments< Dim, P > &a_CA, const IndexedMoments< Dim, P > &a_CB)
Definition: IndexedMomentsImplem.H:418
static void setSize()
Definition: IndexedMomentsImplem.H:383
T sum() const
Definition: IndexTMI.H:260
void checkMoments(IndexedMoments< Dim, ORDER > &a_moments, Vector< IndexTM< int, Dim > > &a_bogusPowers, const Real &a_dx, const Real &a_tolerance, const bool &a_ebMoment, const bool &a_bindMoments)
Definition: IndexedMomentsImplem.H:40
virtual bool ok()
Definition: MomentIterator.H:98
IndexedMoments< Dim, P > & operator+=(const IndexedMoments< Dim, P > &increment)
add compoenentwise
Definition: IndexedMomentsImplem.H:210
Real factorial(const int n)
Calculates factorial for an integer.
Definition: Factorial.H:25
static void setMultiIndicies()
Definition: IndexedMomentsImplem.H:397
size_t size() const
Definition: Vector.H:192
Real getMoment(const IndexTM< int, Dim > &a_mono, const map< IndexTM< int, Dim >, Real > &a_mapin) const
Definition: IndexedMomentsImplem.H:167
void divideByFactorial()
divides each entry by p!
Definition: IndexedMomentsImplem.H:446
#define CH_IM_MAX_POWER
Definition: IndexedMoments.H:23
static void Abort(const char *const a_msg=m_nullString)
Print out message to cerr and exit via abort() (if serial) or MPI_Abort() (if parallel).