Chombo + EB  3.2
MinimalCCCMImplem.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 _MINIMALCCCMIMPLEM_H_
12 #define _MINIMALCCCMIMPLEM_H_
13 
14 #include <cmath>
15 #include <cstdio>
16 #include <cstdlib>
17 #include <fstream>
18 #include <iostream>
19 #include <iomanip>
20 #include <string>
21 
22 #include "MayDay.H"
23 
24 #include "NoRefinement.H"
25 #include "LSProblem.H"
26 #include "Factorial.H"
27 #include "NamespaceHeader.H"
28 
29 // Null constructor
30 template <int dim> MinimalCCCM<dim>::MinimalCCCM()
31 {
32 }
33 
34 // Copy constructor
35 template <int dim> MinimalCCCM<dim>::MinimalCCCM(const MinimalCCCM<dim>& a_MinimalCCCM)
36  :m_cutCellMoments(a_MinimalCCCM.m_cutCellMoments)
37 {
38 }
39 
40 // Use this for initializing
41 template <int dim> MinimalCCCM<dim>::MinimalCCCM(const IFData<dim>& a_info)
42  :m_cutCellMoments(a_info)
43 {
44  m_cutCellMoments.m_bdCCOn = false;
45 
46  for (int hilo = 0; hilo < 2; ++hilo)
47  {
48  for (int idir = 0; idir < dim; ++idir)
49  {
50  // Identifier for which boundary cutCellMoment
51  Iv2 bdId;
52  bdId[BDID_DIR] = idir;
53  bdId[BDID_HILO] = hilo;
54 
55  IFData<dim-1> reducedInfo(a_info,a_info.m_maxOrder+1,idir,hilo);
56 
57  CutCellMoments<dim-1>bdCutCellMoments(reducedInfo);
58 
59  m_cutCellMoments.m_bdCutCellMoments[bdId] = bdCutCellMoments;
60 
61 
62  // Notice whether at least one lower dimensional cutCell is on the
63  // interface
64  if (reducedInfo.m_allVerticesOn)
65  {
66  m_cutCellMoments.m_bdCCOn = true;
67  }
68  }
69  }
70 }
71 
72 // Destructor
73 template <int dim> MinimalCCCM<dim>::~MinimalCCCM()
74 {
75 }
76 
77 
78 // Solve
79 template <int dim> void MinimalCCCM<dim>::computeMoments(const int & a_orderPmax,
80  const int & a_degreePmax)
81 {
82  CH_TIME("computemoments");
83  int zerothOrderOfAccuracy = 0;
84  int highestDegree = a_orderPmax + a_degreePmax;
85 
86  Vector<Real> RNorm(3);
87  for (int i = 0; i < 3; i++)
88  {
89  RNorm[i] = LARGEREALVAL;
90  }
91 
92  for (int i = 0; i <= highestDegree; i++)
93  {
94  m_cutCellMoments.m_residual.push_back(RNorm);
95  }
96 
97  m_boundaryMomentsComputed = false;
98 
99  computeMomentsRecursively(zerothOrderOfAccuracy,
100  highestDegree);
101 
102  IndexTM<int,dim> refineInDir;
103 
104  // Move moments from local coord to parent coord
105  IndexTM<Real,dim> delta = m_cutCellMoments.m_IFData.m_parentCoord.m_origin;
106  delta -= m_cutCellMoments.m_IFData.m_localCoord .m_origin;
107 
108  PthMoment copyMoments = m_cutCellMoments.m_moments;
109  for (typename PthMoment::const_iterator it = copyMoments.begin();
110  it != copyMoments.end(); ++it)
111  {
112  IvDim mono = it->first;
113  m_cutCellMoments.m_moments[mono] = m_cutCellMoments.changeMomentCoordinates(copyMoments, mono, delta);
114  }
115 
116  PthMoment copyEBMoments = m_cutCellMoments.m_EBmoments;
117  for (typename PthMoment::const_iterator it = copyEBMoments.begin(); it != copyEBMoments.end(); ++it)
118  {
119  IvDim mono = it->first;
120  m_cutCellMoments.m_EBmoments[mono] = m_cutCellMoments.changeMomentCoordinates(copyEBMoments, mono, delta);
121  }
122 
123  // Move bdCutCell moments from parent coord to cell center coord. From here on bdCutCell moments will
124  // not be used in any least squares problem
125 
126  for (typename BdCutCellMoments::iterator it = m_cutCellMoments.m_bdCutCellMoments.begin(); it != m_cutCellMoments.m_bdCutCellMoments.end(); ++it)
127  {
128  it->second.changeMomentCoordinatesToCellCenter();
129  }
130 }
131 
132 template <int dim> void MinimalCCCM<dim>::computeMomentsRecursively(const int & a_orderPmax,
133  const int & a_degreePmax)
134 {
135  CH_TIME("computemomentsRecursively");
136  //pout() << "a_orderPmax = " << a_orderPmax << endl;
137  //pout() << "a_degreePmax = " << a_degreePmax << endl;
138  //pout() << "m_cutCellMoments.m_IFData.m_maxOrder " << m_cutCellMoments.m_IFData.m_maxOrder << endl;
139  CH_assert(m_cutCellMoments.m_IFData.m_maxOrder >= a_orderPmax);
140 
141  if (m_cutCellMoments.m_IFData.m_allVerticesOut)
142  {
143  for (int iOrder = 0; iOrder <= a_orderPmax; ++iOrder)
144  {
145  LSProblem<dim> lsp(iOrder + a_degreePmax, false);
146 
147  // Fill moments
148  const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
149  for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
150  {
151  m_cutCellMoments.m_EBmoments[it->first] = 0.0;
152  }
153 
154  const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
155  for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin(); it != monDegreePLess1.end(); ++it)
156  {
157  m_cutCellMoments.m_moments[it->first] = 0.0;
158  }
159  }
160  }
161  else if (m_cutCellMoments.m_IFData.m_allVerticesIn && !m_cutCellMoments.m_bdCCOn)
162  {
163  for (int iOrder = 0; iOrder <= a_orderPmax; ++iOrder)
164  {
165  LSProblem<dim> lsp(iOrder + a_degreePmax, false);
166 
167  // Fill moments of degree P and P-1
168  const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
169  for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
170  {
171  m_cutCellMoments.m_EBmoments[it->first] = 0.0;
172  }
173 
174  const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
175  for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin(); it != monDegreePLess1.end(); ++it)
176  {
177  m_cutCellMoments.m_moments[it->first] = m_cutCellMoments.fullCellQuadrature(it->first,m_cutCellMoments.m_IFData.m_parentCoord);
178  }
179  }
180  }
181  else
182  {
183  // Only compute the boundary moments if they haven't already been
184  // computed (earlier in the recursion)
185  if (!m_boundaryMomentsComputed)
186  {
187  for (typename BdCutCellMoments::iterator it = m_cutCellMoments.m_bdCutCellMoments.begin(); it != m_cutCellMoments.m_bdCutCellMoments.end(); ++it)
188  {
189  MinimalCCCM<dim-1> subProblem(it->second.m_IFData);
190  subProblem.computeMoments(a_orderPmax,a_degreePmax+1);
191  it->second = subProblem.m_cutCellMoments;
192  }
193 
194  m_boundaryMomentsComputed = true;
195  }
196 
197  // Make a LSProb
198  IvDim zeroDerivative = IndexTM<int,dim>::Zero;
199  LSProblem<dim> lsp(a_orderPmax,a_degreePmax, false,m_cutCellMoments.m_IFData.m_normalDerivatives[zeroDerivative]);
200 
201  Vector<Real> rhs = computeRhs(lsp,a_orderPmax);
202 
203  // Solve the problem and return residual
204  int lsCode=lsp.invertNormalEq(rhs,m_cutCellMoments.m_residual[a_degreePmax]);
205  if (lsCode != 0)
206  {
207  pout() << "Geometry generation least squares problem failed with residual:"
208  << m_cutCellMoments.m_residual[a_degreePmax]<< endl;
209  lsp.print(pout());
210  pout () << "Problem occurred generating geometry for these cut cell moments: " << m_cutCellMoments << endl;
211  MayDay::Error("Geometry generation error.[2]");
212  }
213 
214 
215  // Fill moments
216  const PthMomentLoc monDegreePLess1 = lsp.getMonomialLocMapDegreePLess1();
217  for (typename PthMomentLoc::const_iterator it = monDegreePLess1.begin();
218  it != monDegreePLess1.end(); ++it)
219  {
220  m_cutCellMoments.m_moments[it->first] = lsp.getUnknown(it->second + lsp.getNumberDegP());
221  }
222 
223  const PthMomentLoc monDegreeP = lsp.getMonomialLocMapDegreeP();
224  for (typename PthMomentLoc::const_iterator it = monDegreeP.begin(); it != monDegreeP.end(); ++it)
225  {
226  m_cutCellMoments.m_EBmoments[it->first] = lsp.getUnknown(it->second);
227  }
228 
229  }
230 
231  if (a_degreePmax > 0)
232  {
233  computeMomentsRecursively(a_orderPmax + 1,
234  a_degreePmax - 1);
235  }
236 }
237 
239  const int & a_orderPmax)
240 {
241  CH_TIME("computeRHS");
242  // Resize rhs
243  int numEq = dim*a_lsp.getNumberDegP();
244  Vector<Real> rhs(numEq);
245 
246  // For each moment iterate thru bd CutCellMoments incrementing same comp
247  // of rhs
248  const LocPthMoment& locMap = a_lsp.getLocMonomialMapDegreeP();
249  for (typename LocPthMoment::const_iterator it = locMap.begin(); it != locMap.end(); ++it)
250  {
251  int jth = it->first;
252  IvDim mono = it->second;
253 
254  int hiSide = 1;
255  int loSide = 0;
256  Iv2 bdId;
257 
258  for (int idir = 0; idir < dim; ++idir)
259  {
260  // Which lower dimensional monomial corresponds (mono,j)
261  IndexTM<int,dim-1> mono1Less;
262  for (int jdir = 0; jdir < dim; ++jdir)
263  {
264  if (jdir < idir)
265  {
266  mono1Less[jdir] = mono[jdir];
267  }
268  else if (jdir > idir)
269  {
270  mono1Less[jdir-1] = mono[jdir];
271  }
272  }
273 
274  bdId[0] = idir;
275  bdId[1] = hiSide;
276 
277  Real hiMom = m_cutCellMoments.m_bdCutCellMoments[bdId].m_moments[mono1Less];
278 
279  bdId[1] = loSide;
280 
281  Real loMom = m_cutCellMoments.m_bdCutCellMoments[bdId].m_moments[mono1Less];
282  int exponent = it->second[idir];
283 
284  Real loSideValue;
285  Real hiSideValue;
286 
287  loSideValue = m_cutCellMoments.m_IFData.m_localCoord.convertDir(-0.5*m_cutCellMoments.m_IFData.m_cellCenterCoord.m_dx[idir],
288  m_cutCellMoments.m_IFData.m_cellCenterCoord,
289  idir);
290 
291  hiSideValue = m_cutCellMoments.m_IFData.m_localCoord.convertDir( 0.5*m_cutCellMoments.m_IFData.m_cellCenterCoord.m_dx[idir],
292  m_cutCellMoments.m_IFData.m_cellCenterCoord,
293  idir);
294 
295  Real loFactor = POW(loSideValue,exponent);
296  Real hiFactor = POW(hiSideValue,exponent);
297 
298  rhs[(dim*jth) + idir] = hiMom*hiFactor - loMom*loFactor;
299 
300  // Add the Taylor series terms
301  for (int order = 1; order <= a_orderPmax; order++)
302  {
303  Vector<IvDim> taylorMonomials;
304 
305  generateMultiIndices(taylorMonomials,order);
306 
307  for (int i = 0; i < taylorMonomials.size(); i++)
308  {
309  const IvDim & taylorMonomial = taylorMonomials[i];
310 
311  IvDim totalMonomial = mono + taylorMonomial;
312 
313  if (m_cutCellMoments.m_EBmoments.find(totalMonomial) !=
314  m_cutCellMoments.m_EBmoments.end())
315  {
316  Real normalDerivative = m_cutCellMoments.m_IFData.m_normalDerivatives[taylorMonomial][idir];
317  Real fact = factorial(taylorMonomial);
318 
319  Real moment = m_cutCellMoments.m_EBmoments[totalMonomial];
320 
321  rhs[(dim*jth) + idir] += normalDerivative * moment / fact;
322  }
323  else
324  {
325  MayDay::Error("Unable to find needed monomial for Taylor series");
326  }
327  }
328  }
329  }
330  }
331 
332  return rhs;
333 }
334 template <int dim> void MinimalCCCM<dim>::print(ostream& a_out) const
335 {
336  m_cutCellMoments.print(a_out);
337 }
338 
339 template <int dim> void MinimalCCCM<dim>::dump() const
340 {
341  print(pout());
342 }
343 
344 // Operators
345 template <int dim> void MinimalCCCM<dim>::operator=(const MinimalCCCM<dim> & a_MinimalCCCM)
346 {
347  // Only copy if the objects are distinct
348  if (this != &a_MinimalCCCM)
349  {
350  m_cutCellMoments = a_MinimalCCCM.m_cutCellMoments;
351  }
352 }
353 
354 template <int dim> ostream& operator<<(ostream & a_out,
355  const MinimalCCCM<dim> & a_MinimalCCCM)
356 {
357  a_MinimalCCCM.print(a_out);
358  return a_out;
359 }
360 
361 template <int dim> Real MinimalCCCM<dim>::factorial(const IvDim & a_multiIndex) const
362 {
363  Real fact = 1.0;
364 
365  for (int i = 0; i < dim; i++)
366  {
367  for (int j = 2; j <= a_multiIndex[i]; j++)
368  {
369  fact *= j;
370  }
371  }
372 
373  return fact;
374 }
375 
376 #include "NamespaceFooter.H"
377 
378 #endif
std::ostream & pout()
Use this in place of std::cout for program output.
map< int, IvDim > LocPthMoment
Definition: MinimalCCCM.H:42
#define BDID_HILO
Definition: Notation.H:89
Real factorial(const IvDim &a_multiIndex) const
Definition: MinimalCCCMImplem.H:361
MinimalCCCM()
Definition: MinimalCCCMImplem.H:30
#define CH_assert(cond)
Definition: CHArray.H:37
#define LARGEREALVAL
Definition: Notation.H:77
CutCellMoments< dim > m_cutCellMoments
Definition: MinimalCCCM.H:78
Definition: MinimalCCCM.H:31
Real getUnknown(int loc)
Definition: LSProblem.H:104
Definition: ComputeCutCellMoments.H:35
Real POW(const Real &a_x, const int &a_p)
computes x^p
Definition: Factorial.H:33
ostream & operator<<(ostream &a_out, const MinimalCCCM< dim > &a_MinimalCCCM)
Definition: MinimalCCCMImplem.H:354
Definition: IndexTM.H:36
map< IvDim, Real > PthMoment
Definition: MinimalCCCM.H:37
Vector< Real > computeRhs(LSProblem< dim > &a_lsp, const int &a_order)
Definition: MinimalCCCMImplem.H:238
void print(ostream &a_out) const
Definition: LSProblemImplem.H:560
void dump() const
Definition: MinimalCCCMImplem.H:339
#define CH_TIME(name)
Definition: CH_Timer.H:82
Definition: IFData.H:42
#define BDID_DIR
Definition: Notation.H:88
void operator=(const MinimalCCCM< dim > &a_MinimalCCCM)
Definition: MinimalCCCMImplem.H:345
double Real
Definition: REAL.H:33
const PthMomentLoc & getMonomialLocMapDegreeP() const
Definition: LSProblem.H:94
int invertNormalEq(const Vector< Real > &a_rhs, Vector< Real > &a_residual)
Definition: LSProblemImplem.H:317
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
void print(ostream &out) const
Definition: MinimalCCCMImplem.H:334
int getNumberDegP()
Definition: LSProblem.H:109
Real factorial(const int n)
Calculates factorial for an integer.
Definition: Factorial.H:25
size_t size() const
Definition: Vector.H:192
void computeMoments(const int &a_orderPmax, const int &a_degreePmax)
Definition: MinimalCCCMImplem.H:79
const PthMomentLoc & getMonomialLocMapDegreePLess1() const
Definition: LSProblem.H:99
int m_maxOrder
Definition: IFData.H:67
~MinimalCCCM()
Definition: MinimalCCCMImplem.H:73
void generateMultiIndices(Vector< IndexTM< int, dim > > &a_indices, const int &a_magnitude)
Definition: MultiIndexImplem.H:18
int dim
Definition: EBInterface.H:146
map< IvDim, int > PthMomentLoc
Definition: MinimalCCCM.H:43
Definition: CutCellMoments.H:32
void computeMomentsRecursively(const int &a_orderPmax, const int &a_degreePmax)
Definition: MinimalCCCMImplem.H:132
const LocPthMoment & getLocMonomialMapDegreeP() const
Definition: LSProblem.H:89