Chombo + EB  3.0
CoefficientInterpolator.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 _COEFFICIENTINTERPOLATOR_H__
12 #define _COEFFICIENTINTERPOLATOR_H__
13 
14 #include "LevelData.H"
15 #include "DataIterator.H"
16 #include "DebugOut.H"
17 #include "BoxIterator.H"
18 #include "NamespaceHeader.H"
19 
20 // Jeffrey Johnson (jnjohnson@lbl.gov)
21 
22 //! \class CoefficientInterpolator
23 //! This abstract base class provides an interface to time-dependent spatial coefficient
24 //! data for various partial differential equations. One obtains constant access to
25 //! coefficient data by specifying the desired time centering. Subclasses of this
26 //! base class define how that coefficient data is computed.
27 //! \tparam LevelDataType The LevelData type used to store coefficient data.
28 template <typename LevelData_, typename SolutionLevelData_ = LevelData_>
30 {
31  public:
32 
33  typedef LevelData_ LevelDataType;
34  typedef SolutionLevelData_ SolutionLevelDataType;
35 
36  //! Base class constructor. Called by every subclass.
37  //! \param a_numComps The number of components in the interpolated coefficient.
38  explicit CoefficientInterpolator(int a_numComps);
39 
40  //! Destructor.
41  virtual ~CoefficientInterpolator();
42 
43  //! Returns the number of components in the interpolated coefficient.
44  int numComps() const;
45 
46  //! Interpolate the coefficient at the given time, placing the result in the given
47  //! LevelData object. This method must be overridden by subclasses.
48  //! \param a_result The LevelData object that will store the result.
49  //! \param a_time The time at which the coefficient is to be evaluated.
50  virtual void interpolate(LevelDataType& a_result,
51  Real a_time);
52 
53  //! Interpolate the coefficient at the given time, placing the result in the given
54  //! LevelData object. This method assumes that the coefficient depends upon the
55  //! solution to the partial differential equation in question, so the solution
56  //! is passed into it as an argument.
57  //! \param a_result The LevelData object that will store the result.
58  //! \param a_solution The solution to the partial differential equation.
59  //! \param a_time The time at which the coefficient is to be evaluated.
60  virtual void interpolate(LevelDataType& a_result,
61  const SolutionLevelDataType& a_solution,
62  Real a_time);
63 
64  //! Returns true if the coefficient depends on the solution to the
65  //! partial differential equation (rendering it nonlinear), false
66  //! otherwise. By default, the coefficient is assumed not to depend
67  //! upon the solution.
68  virtual bool dependsUponSolution() const;
69 
70  //! Computes the derivative of the coefficient with respect to the
71  //! solution at the desired time. By default, this sets \a a_deriv to 0.
72  //! \param a_prime The coefficient derivative data will be stored here.
73  //! \param a_solution The solution to the partial differential equation.
74  //! \param a_time The time at which to compute the coefficient data.
75  virtual void interpolatePrime(LevelDataType& a_prime,
76  const SolutionLevelDataType& a_solution,
77  Real a_time);
78 
79  //! This virtual void method performs the iterative nonlinear solve
80  //! \f$A(\phi) \phi - f(\vec{x}) = 0\f$ for \f$\phi\f$.
81  //! \param a_phi The solution to the equation, \f$\phi\f$, will be stored here.
82  //! \param a_f The term \f$f(\vec{x})\f$ in the equation.
83  //! \param a_time The time at which the equation is solved.
84  //! \param a_phi0 The initial estimate for \f$\phi\f$.
85  //! \param a_tolerance The threshold for the error in the equation that
86  //! dictates when iteration should cease.
87  virtual void solve(SolutionLevelDataType& a_phi,
88  const SolutionLevelDataType& a_f,
89  Real a_time,
90  const SolutionLevelDataType& a_phi0,
91  Real a_tolerance);
92 
93  //! This helper method solves the nonlinear equation
94  //! \f$A(\phi) \phi - f(\vec{x}) = 0\f$ for \f$\phi\f$ using Newton-Raphson
95  //! iteration. It can be used by subclasses to implement the solve()
96  //! virtual method.
97  //! \param a_phi The solution to the equation, \f$\phi\f$, will be stored here.
98  //! \param a_f The term \f$f(\vec{x})\f$ in the equation.
99  //! \param a_time The time at which the equation is solved.
100  //! \param a_phi0 The initial estimate for \f$\phi\f$.
101  //! \param a_tolerance The threshold for the error in the equation that
102  //! dictates when iteration should cease.
104  const SolutionLevelDataType& a_f,
105  Real a_time,
106  const SolutionLevelDataType& a_phi0,
107  Real a_tolerance);
108 
109  private:
110 
111  // Number of components.
113 
114  // Prevents tail-chasing.
115  bool m_inCall;
116 
117  // Disallowed operators.
121 };
122 
123 
124 //-----------------------------------------------------------------------
125 template <typename LevelData_, typename SolutionLevelData_>
127 CoefficientInterpolator(int a_numComps)
128  :m_numComps(a_numComps),
129  m_inCall(false)
130 {
131 }
132 //-----------------------------------------------------------------------
133 
134 //-----------------------------------------------------------------------
135 template <typename LevelData_, typename SolutionLevelData_>
138 {
139 }
140 //-----------------------------------------------------------------------
141 
142 //-----------------------------------------------------------------------
143 template <typename LevelData_, typename SolutionLevelData_>
144 int
146 numComps() const
147 {
148  return m_numComps;
149 }
150 //-----------------------------------------------------------------------
151 
152 //-----------------------------------------------------------------------
153 template <typename LevelData_, typename SolutionLevelData_>
154 bool
157 {
158  return false;
159 }
160 //-----------------------------------------------------------------------
161 
162 //-----------------------------------------------------------------------
163 template <typename LevelData_, typename SolutionLevelData_>
164 void
166 interpolate(LevelData_& a_result,
167  Real a_time)
168 {
169  if (dependsUponSolution())
170  {
171  MayDay::Error("The solution must be passed to this interpolator!");
172  }
173  if (m_inCall)
174  {
175  MayDay::Error("Neither the solution-independent nor -dependent interpolate method has\n"
176  "been defined for this subclass.");
177  }
178 
179  // It doesn't matter what we pass in for the solution.
180  m_inCall = true;
181  SolutionLevelData_ phoneyBaloney(a_result.disjointBoxLayout(), 1, a_result.ghostVect());
182  interpolate(a_result, phoneyBaloney, a_time);
183  m_inCall = false;
184 }
185 //-----------------------------------------------------------------------
186 
187 //-----------------------------------------------------------------------
188 template <typename LevelData_, typename SolutionLevelData_>
189 void
191 interpolate(LevelData_& a_result,
192  const SolutionLevelData_& a_solution,
193  Real a_time)
194 {
195  if (m_inCall)
196  {
197  MayDay::Error("Neither the solution-independent nor -dependent interpolate method has\n"
198  "been defined for this subclass.");
199  }
200  if (!dependsUponSolution())
201  {
202  m_inCall = true;
203  interpolate(a_result, a_time);
204  m_inCall = false;
205  }
206  else
207  {
208  MayDay::Error("The solution-dependent interpolate method has not been defined for this subclass.");
209  }
210 }
211 //-----------------------------------------------------------------------
212 
213 //-----------------------------------------------------------------------
214 template <typename LevelData_, typename SolutionLevelData_>
215 void
217 interpolatePrime(LevelData_& a_prime,
218  const SolutionLevelData_& a_solution,
219  Real a_time)
220 {
221  if (!dependsUponSolution())
222  {
223  for (DataIterator dit = a_prime.disjointBoxLayout().dataIterator(); dit.ok(); ++dit)
224  a_prime[dit()].setVal(0.0);
225  }
226  else
227  MayDay::Error("The derivative w.r.t. the solution has not been implemented for this subclass.");
228 }
229 //-----------------------------------------------------------------------
230 
231 //-----------------------------------------------------------------------
232 template <typename LevelData_, typename SolutionLevelData_>
233 void
235 solve(SolutionLevelData_& a_phi,
236  const SolutionLevelData_& a_f,
237  Real a_time,
238  const SolutionLevelData_& a_phi0,
239  Real a_tolerance)
240 {
241  MayDay::Error("CoefficientInterpolator::solve() must be implemented for this subclass!");
242 }
243 //-----------------------------------------------------------------------
244 
245 //-----------------------------------------------------------------------
246 template <typename LevelData_, typename SolutionLevelData_>
247 void
249 NewtonRaphson(SolutionLevelData_& a_phi,
250  const SolutionLevelData_& a_f,
251  Real a_time,
252  const SolutionLevelData_& a_phi0,
253  Real a_tolerance)
254 {
255  // How bad is the initial estimate?
256  Real maxError = -FLT_MAX;
257  LevelData_ F(a_phi0.disjointBoxLayout(), a_phi0.nComp()),
258  A(a_phi0.disjointBoxLayout(), a_phi0.nComp());
259  interpolate(A, a_phi0, a_time);
260  for (DataIterator dit = a_phi0.dataIterator(); dit.ok(); ++dit)
261  {
262  a_phi[dit()].copy(a_phi0[dit()]); // phi <- phi0.
263  F[dit()].copy(A[dit()]);
264  F[dit()] *= a_phi[dit()];
265  F[dit()] -= a_f[dit()];
266  maxError = Max(Abs(F[dit()].max()), Abs(F[dit()].min()));
267  }
268  if (maxError < a_tolerance)
269  return; // Actually, our initial estimate was fine!
270 
271  LevelData_ Aprime(a_phi.disjointBoxLayout(), a_phi.nComp());
272  while (maxError > a_tolerance)
273  {
274  // Compute the correction term F/F', apply it to phi, and recompute F.
275  // FIXME: This is slow and probably should be moved to Fortran.
276  interpolatePrime(Aprime, a_phi, a_time);
277  for (DataIterator dit = a_phi.dataIterator(); dit.ok(); ++dit)
278  {
279  Box box = a_phi.disjointBoxLayout().get(dit());
280  const FArrayBox& a = A[dit()];
281  const FArrayBox& aprime = Aprime[dit()];
282  const FArrayBox& f = F[dit()];
283  FArrayBox& phi = a_phi[dit()];
284  for (BoxIterator bit(box); bit.ok(); ++bit)
285  {
286  IntVect i = bit();
287  for (int n = 0; n < F.nComp(); ++n)
288  {
289  if (Abs(f(i, n)) > a_tolerance)
290  {
291  Real fprime = aprime(i, n) * phi(i, n) + a(i, n);
292  phi(i, n) -= f(i, n) / fprime;
293  }
294  }
295  }
296  }
297 
298  // Recompute F with the corrected value of phi.
299  interpolate(A, a_phi, a_time);
300  for (DataIterator dit = a_phi.dataIterator(); dit.ok(); ++dit)
301  {
302  // Recompute F.
303  F[dit()].copy(A[dit()]);
304  F[dit()] *= a_phi[dit()];
305  F[dit()] -= a_f[dit()];
306 
307  // Measure the error again.
308  maxError = Max(Abs(F[dit()].max()), Abs(F[dit()].min()));
309  }
310  }
311 }
312 //-----------------------------------------------------------------------
313 
314 #include "NamespaceFooter.H"
315 
316 #endif
bool ok()
Definition: BoxIterator.H:215
LevelData_ LevelDataType
Definition: CoefficientInterpolator.H:33
IndexTM< T, N > min(const IndexTM< T, N > &a_p1, const IndexTM< T, N > &a_p2)
Definition: IndexTMI.H:396
int numComps() const
Returns the number of components in the interpolated coefficient.
Definition: CoefficientInterpolator.H:146
Definition: DataIterator.H:140
iterates through the IntVects of a Box
Definition: BoxIterator.H:37
Definition: EBInterface.H:45
Definition: CoefficientInterpolator.H:29
double Real
Definition: REAL.H:33
virtual void interpolate(LevelDataType &a_result, Real a_time)
Definition: CoefficientInterpolator.H:166
int m_numComps
Definition: CoefficientInterpolator.H:112
IndexTM< T, N > max(const IndexTM< T, N > &a_p1, const IndexTM< T, N > &a_p2)
Definition: IndexTMI.H:403
T Abs(const T &a_a)
Definition: Misc.H:53
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.
bool m_inCall
Definition: CoefficientInterpolator.H:115
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
SolutionLevelData_ SolutionLevelDataType
Definition: CoefficientInterpolator.H:34
virtual ~CoefficientInterpolator()
Destructor.
Definition: CoefficientInterpolator.H:137
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
Definition: FArrayBox.H:44
virtual void solve(SolutionLevelDataType &a_phi, const SolutionLevelDataType &a_f, Real a_time, const SolutionLevelDataType &a_phi0, Real a_tolerance)
Definition: CoefficientInterpolator.H:235
T Max(const T &a_a, const T &a_b)
Definition: Misc.H:39
virtual bool ok() const
return true if this iterator is still in its Layout
Definition: LayoutIterator.H:110
virtual bool dependsUponSolution() const
Definition: CoefficientInterpolator.H:156
virtual void interpolatePrime(LevelDataType &a_prime, const SolutionLevelDataType &a_solution, Real a_time)
Definition: CoefficientInterpolator.H:217
void NewtonRaphson(SolutionLevelDataType &a_phi, const SolutionLevelDataType &a_f, Real a_time, const SolutionLevelDataType &a_phi0, Real a_tolerance)
Definition: CoefficientInterpolator.H:249
CoefficientInterpolator & operator=(const CoefficientInterpolator &)