Chombo + EB + MF  3.2
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 ([email protected])
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.
103  void NewtonRaphson(SolutionLevelDataType& a_phi,
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 template <typename LevelData_, typename SolutionLevelData_>
126 CoefficientInterpolator(int a_numComps)
127  :m_numComps(a_numComps),
128  m_inCall(false)
129 {
130 }
131 //-----------------------------------------------------------------------
132 
133 //-----------------------------------------------------------------------
134 template <typename LevelData_, typename SolutionLevelData_>
137 {
138 }
139 //-----------------------------------------------------------------------
140 
141 //-----------------------------------------------------------------------
142 template <typename LevelData_, typename SolutionLevelData_>
143 int
145 numComps() const
146 {
147  return m_numComps;
148 }
149 //-----------------------------------------------------------------------
150 
151 //-----------------------------------------------------------------------
152 template <typename LevelData_, typename SolutionLevelData_>
153 bool
156 {
157  return false;
158 }
159 //-----------------------------------------------------------------------
160 
161 //-----------------------------------------------------------------------
162 template <typename LevelData_, typename SolutionLevelData_>
163 void
165 interpolate(LevelData_& a_result,
166  Real a_time)
167 {
168  if (dependsUponSolution())
169  {
170  MayDay::Error("The solution must be passed to this interpolator!");
171  }
172  if (m_inCall)
173  {
174  MayDay::Error("Neither the solution-independent nor -dependent interpolate method has\n"
175  "been defined for this subclass.");
176  }
177 
178  // It doesn't matter what we pass in for the solution.
179  m_inCall = true;
180  SolutionLevelData_ phoneyBaloney(a_result.disjointBoxLayout(), 1, a_result.ghostVect());
181  interpolate(a_result, phoneyBaloney, a_time);
182  m_inCall = false;
183 }
184 //-----------------------------------------------------------------------
185 
186 //-----------------------------------------------------------------------
187 template <typename LevelData_, typename SolutionLevelData_>
188 void
190 interpolate(LevelData_& a_result,
191  const SolutionLevelData_& a_solution,
192  Real a_time)
193 {
194  if (m_inCall)
195  {
196  MayDay::Error("Neither the solution-independent nor -dependent interpolate method has\n"
197  "been defined for this subclass.");
198  }
199  if (!dependsUponSolution())
200  {
201  m_inCall = true;
202  interpolate(a_result, a_time);
203  m_inCall = false;
204  }
205  else
206  {
207  MayDay::Error("The solution-dependent interpolate method has not been defined for this subclass.");
208  }
209 }
210 //-----------------------------------------------------------------------
211 
212 //-----------------------------------------------------------------------
213 template <typename LevelData_, typename SolutionLevelData_>
214 void
216 interpolatePrime(LevelData_& a_prime,
217  const SolutionLevelData_& a_solution,
218  Real a_time)
219 {
220  if (!dependsUponSolution())
221  {
222  for (DataIterator dit = a_prime.disjointBoxLayout().dataIterator(); dit.ok(); ++dit)
223  a_prime[dit()].setVal(0.0);
224  }
225  else
226  MayDay::Error("The derivative w.r.t. the solution has not been implemented for this subclass.");
227 }
228 //-----------------------------------------------------------------------
229 
230 //-----------------------------------------------------------------------
231 template <typename LevelData_, typename SolutionLevelData_>
232 void
234 solve(SolutionLevelData_& a_phi,
235  const SolutionLevelData_& a_f,
236  Real a_time,
237  const SolutionLevelData_& a_phi0,
238  Real a_tolerance)
239 {
240  MayDay::Error("CoefficientInterpolator::solve() must be implemented for this subclass!");
241 }
242 //-----------------------------------------------------------------------
243 
244 //-----------------------------------------------------------------------
245 template <typename LevelData_, typename SolutionLevelData_>
246 void
248 NewtonRaphson(SolutionLevelData_& a_phi,
249  const SolutionLevelData_& a_f,
250  Real a_time,
251  const SolutionLevelData_& a_phi0,
252  Real a_tolerance)
253 {
254  // How bad is the initial estimate?
255  Real maxError = -FLT_MAX;
256  LevelData_ F(a_phi0.disjointBoxLayout(), a_phi0.nComp()),
257  A(a_phi0.disjointBoxLayout(), a_phi0.nComp());
258  interpolate(A, a_phi0, a_time);
259  for (DataIterator dit = a_phi0.dataIterator(); dit.ok(); ++dit)
260  {
261  a_phi[dit()].copy(a_phi0[dit()]); // phi <- phi0.
262  F[dit()].copy(A[dit()]);
263  F[dit()] *= a_phi[dit()];
264  F[dit()] -= a_f[dit()];
265  maxError = Max(Abs(F[dit()].max()), Abs(F[dit()].min()));
266  }
267  if (maxError < a_tolerance)
268  return; // Actually, our initial estimate was fine!
269 
270  LevelData_ Aprime(a_phi.disjointBoxLayout(), a_phi.nComp());
271  while (maxError > a_tolerance)
272  {
273  // Compute the correction term F/F', apply it to phi, and recompute F.
274  // FIXME: This is slow and probably should be moved to Fortran.
275  interpolatePrime(Aprime, a_phi, a_time);
276  for (DataIterator dit = a_phi.dataIterator(); dit.ok(); ++dit)
277  {
278  Box box = a_phi.disjointBoxLayout().get(dit());
279  const FArrayBox& a = A[dit()];
280  const FArrayBox& aprime = Aprime[dit()];
281  const FArrayBox& f = F[dit()];
282  FArrayBox& phi = a_phi[dit()];
283  for (BoxIterator bit(box); bit.ok(); ++bit)
284  {
285  IntVect i = bit();
286  for (int n = 0; n < F.nComp(); ++n)
287  {
288  if (Abs(f(i, n)) > a_tolerance)
289  {
290  Real fprime = aprime(i, n) * phi(i, n) + a(i, n);
291  phi(i, n) -= f(i, n) / fprime;
292  }
293  }
294  }
295  }
296 
297  // Recompute F with the corrected value of phi.
298  interpolate(A, a_phi, a_time);
299  for (DataIterator dit = a_phi.dataIterator(); dit.ok(); ++dit)
300  {
301  // Recompute F.
302  F[dit()].copy(A[dit()]);
303  F[dit()] *= a_phi[dit()];
304  F[dit()] -= a_f[dit()];
305 
306  // Measure the error again.
307  maxError = Max(Abs(F[dit()].max()), Abs(F[dit()].min()));
308  }
309  }
310 }
311 //-----------------------------------------------------------------------
312 
313 #include "NamespaceFooter.H"
314 
315 #endif
bool ok()
Definition: BoxIterator.H:281
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:394
virtual bool ok() const
return true if this iterator is still in its Layout
Definition: LayoutIterator.H:117
Definition: DataIterator.H:190
iterates through the IntVects of a Box
Definition: BoxIterator.H:37
Definition: EBInterface.H:45
void const char const int const int const int const Real const Real * A
Definition: Lapack.H:83
Definition: CoefficientInterpolator.H:29
double Real
Definition: REAL.H:33
virtual void interpolate(LevelDataType &a_result, Real a_time)
Definition: CoefficientInterpolator.H:165
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:401
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.
virtual bool dependsUponSolution() const
Definition: CoefficientInterpolator.H:155
int numComps() const
Returns the number of components in the interpolated coefficient.
Definition: CoefficientInterpolator.H:145
bool m_inCall
Definition: CoefficientInterpolator.H:115
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:469
SolutionLevelData_ SolutionLevelDataType
Definition: CoefficientInterpolator.H:34
virtual ~CoefficientInterpolator()
Destructor.
Definition: CoefficientInterpolator.H:136
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
Definition: FArrayBox.H:45
virtual void solve(SolutionLevelDataType &a_phi, const SolutionLevelDataType &a_f, Real a_time, const SolutionLevelDataType &a_phi0, Real a_tolerance)
Definition: CoefficientInterpolator.H:234
T Max(const T &a_a, const T &a_b)
Definition: Misc.H:39
virtual void interpolatePrime(LevelDataType &a_prime, const SolutionLevelDataType &a_solution, Real a_time)
Definition: CoefficientInterpolator.H:216
void NewtonRaphson(SolutionLevelDataType &a_phi, const SolutionLevelDataType &a_f, Real a_time, const SolutionLevelDataType &a_phi0, Real a_tolerance)
Definition: CoefficientInterpolator.H:248
CoefficientInterpolator & operator=(const CoefficientInterpolator &)