LifeV
analyticalSol.hpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 #ifndef ANALYTICALSOL_H_
28 #define ANALYTICALSOL_H_
29 
30 #include <lifev/core/LifeV.hpp>
31 
32 namespace LifeV
33 {
34 
35 /*!
36  \typedef The differential problem at hand
37 */
38 enum ADRProblemSolution
39 {
40  POISSON_POLYNOMIAL, //!< A polynomial solution to the Poisson problem
41  POISSON_TRIGONOMETRIC, //!< A trigonometric solution to the Poisson problem
42  ADR_STEADY_POLYNOMIAL, //!< A polynomial solution to a general ADR problem
43  ADR_UNSTEADY_POLYNOMIAL //!< A polynomial solution to a general ADR problem
44 };
45 
46 
47 // Coefficient for Robin boundary condition
48 Real coefRobin (const Real& /*t*/, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const ID& /*i*/)
49 {
50  return 2.0;
51 }
52 
53 #ifdef TWODIM
54 
55 class AnalyticalSol
56 {
57 public:
58  // Give to the class a "functor" interface
59  Real operator() (const Real& t, const Real& x, const Real& y, const Real& z, const UInt& ic) const
60  {
61  return u_ex ( t, x, y, z, ic );
62  }
63  // This method is required by the error calculation routines
64  Real grad (UInt icoor, const Real& t, const Real& x, const Real& y, const Real& z, const UInt& ic) const
65  {
66  return grad_ex ( icoor, t, x, y, z, ic );
67  }
68  // The exact solution to the problem
69  static Real u_ex (const Real& /*t*/, const Real& x, const Real& y, const Real& /*z*/, const UInt& /*ic=0*/)
70  {
71  // return x*x + y*y; // Polynomial solution
72  return std::sin (2 * Pi * x) * std::sin (2 * Pi * y); // Trigonometric solution
73  }
74  // The gradient of the exact solution
75  static Real grad_ex (const UInt& icoor, const Real& /*t*/, const Real& x, const Real& y, const Real& /*z*/, const UInt& /*ic=0*/)
76  {
77  switch (icoor)
78  {
79  case 1: //der_x
80  // return 2*x; // Polynomial solution
81  return 2 * Pi * std::cos (2 * Pi * x) * std::sin (2 * Pi * y); // Trigonometric solution
82  case 2: //der_y
83  // return 2*x; // Polynomial solution
84  return 2 * Pi * std::sin (2 * Pi * x) * std::cos (2 * Pi * y); // Trigonometric solution
85  default:
86  return 0;
87  }
88  }
89  // The source term
90  static Real source (const Real& /*t*/, const Real& x, const Real& y, const Real& /*z*/, const UInt& /*ic=0*/)
91  {
92  // return -4; // Polynomial solution
93  return 8 * Pi * Pi * std::sin (2 * Pi * x) * std::sin (2 * Pi * y); // Trigonometric solution
94  }
95  // The normal derivative to the domain boundary
96  // BEWARE: this is case dependent! (works for cubic domains)
97  static Real fNeumann (const Real& t, const Real& x, const Real& y, const Real& z, const UInt& ic)
98  {
99  Real n[2] = {0, 0};
100  Real out = 0;
101  if ( x == xMin )
102  {
103  n[0] = -1.;
104  }
105  else if ( x == xMax )
106  {
107  n[0] = 1.;
108  }
109  else if ( y == yMin )
110  {
111  n[1] = -1.;
112  }
113  else if ( y == yMax )
114  {
115  n[1] = 1.;
116  }
117  else
118  {
119  std::cout << "strange point: x=" << x << " y=" << y << " z=" << z
120  << std::endl;
121  }
122 
123  for (UInt k = 0; k < 2; k++) //mu gradu n
124  {
125  out += grad_ex (k + 1, t, x, y, z, ic) * n[k];
126  }
127 
128  return out;
129  }
130  // Robin boundary condition compatible with the exact solution
131  static Real fRobin (const Real& t, const Real& x, const Real& y, const Real& z, const UInt& ic)
132  {
133  return fNeumann (t, x, y, z, ic) + coefRobin (t, x, y, z, ic) * u_ex (t, x, y, z, ic);
134  }
135  // Initialization of the private members
136  static void setup ( const GetPot& dataFile, const std::string& section = "poisson" )
137  {
138  // We want a slash dividing the data file section from the variable name but only
139  // when not provided by the user or when not looking in the root of the data file
140  std::string corrected_section ( section );
141  if ( ( ! section.empty() ) && ( section[section.length() - 1] != '/' ) )
142  {
143  corrected_section = section + '/';
144  }
145 
146  xMin = dataFile ( (corrected_section + "problem/xMin").data(), -1. ) ;
147  xMax = dataFile ( (corrected_section + "problem/xMax").data(), 1. ) ;
148  yMin = dataFile ( (corrected_section + "problem/yMin").data(), -1. ) ;
149  yMax = dataFile ( (corrected_section + "problem/yMax").data(), 1. ) ;
150  }
151 
152 private:
153  // The dimensions of the domain
154  static Real xMin;
155  static Real xMax;
156  static Real yMin;
157  static Real yMax;
158 
159 };
160 
161 // Declaration and default initialization of the static variables
162 Real AnalyticalSol::xMin (0.);
163 Real AnalyticalSol::xMax (1.);
164 Real AnalyticalSol::yMin (0.);
165 Real AnalyticalSol::yMax (1.);
166 
167 
168 
169 Real beta ( const Real& /* t */,
170  const Real& ,
171  const Real& ,
172  const Real& ,
173  const ID& icomp )
174 {
175  switch (icomp)
176  {
177  case 1:
178  return 0;
179  case 2:
180  return 0;
181  default:
182  return 0;
183  }
184 }
185 
186 
187 #elif defined THREEDIM
188 
189 class AnalyticalSol
190 {
191 public:
192  // Give to the class a "functor" interface
193  Real operator() (const Real& t, const Real& x, const Real& y, const Real& z, const UInt& ic) const
194  {
195  return u_ex ( t, x, y, z, ic );
196  }
197  // This method is required by the error calculation routines
198  Real grad (UInt icoor, const Real& t, const Real& x, const Real& y, const Real& z, const UInt& ic) const
199  {
200  return grad_ex ( icoor, t, x, y, z, ic );
201  }
202  // The exact solution to the problem
203  static Real u_ex (const Real& t, const Real& x, const Real& y, const Real& z, const UInt& /*ic*/)
204  {
205  switch ( problemSolution )
206  {
207  case POISSON_POLYNOMIAL: // Poisson problem, polynomial solution
208  case ADR_STEADY_POLYNOMIAL: // ADR problem, polynomial solution
209  return x * x + y * y + z * z;
210  case POISSON_TRIGONOMETRIC: // Poisson problem, trigonometric solution
211  return std::sin (2 * Pi * x) * std::sin (2 * Pi * y) * std::sin (2 * Pi * z);
212  case ADR_UNSTEADY_POLYNOMIAL: // ADR problem, polynomial solution
213  return t * t + x * x + y * y + z * z;
214  default:
215  return 0;
216  }
217  }
218  // The gradient of the exact solution
219  static Real grad_ex (const UInt& icoor,
220  const Real& /*t*/, const Real& x, const Real& y, const Real& z, const UInt& /*ic*/)
221  {
222  switch (icoor)
223  {
224  case 1: //der_x
225  switch ( problemSolution )
226  {
227  case POISSON_POLYNOMIAL: // Poisson problem, polynomial solution
228  case ADR_STEADY_POLYNOMIAL: // ADR problem, polynomial solution
229  case ADR_UNSTEADY_POLYNOMIAL: // ADR problem, polynomial solution
230  return 2 * x; // Polynomial solution
231  case POISSON_TRIGONOMETRIC: // Poisson problem, trigonometric solution
232  return 2 * Pi * std::cos (2 * Pi * x) * std::sin (2 * Pi * y) * std::sin (2 * Pi * z); // Trigonometric solution
233  }
234  case 2: //der_y
235  switch ( problemSolution )
236  {
237  case POISSON_POLYNOMIAL: // Poisson problem, polynomial solution
238  case ADR_STEADY_POLYNOMIAL: // ADR problem, polynomial solution
239  case ADR_UNSTEADY_POLYNOMIAL: // ADR problem, polynomial solution
240  return 2 * y; // Polynomial solution
241  case POISSON_TRIGONOMETRIC: // Poisson problem, trigonometric solution
242  return 2 * Pi * std::sin (2 * Pi * x) * std::cos (2 * Pi * y) * std::sin (2 * Pi * z); // Trigonometric solution
243  }
244  case 3: //der_z
245  switch ( problemSolution )
246  {
247  case POISSON_POLYNOMIAL: // Poisson problem, polynomial solution
248  case ADR_STEADY_POLYNOMIAL: // ADR problem, polynomial solution
249  case ADR_UNSTEADY_POLYNOMIAL: // ADR problem, polynomial solution
250  return 2 * z; // Polynomial solution
251  case POISSON_TRIGONOMETRIC: // Poisson problem, trigonometric solution
252  return 2 * Pi * std::sin (2 * Pi * x) * std::sin (2 * Pi * y) * std::cos (2 * Pi * z); // Trigonometric solution
253  }
254  default:
255  return 0;
256  }
257  }
258  // The source term
259  static Real fSource (const Real& t, const Real& x, const Real& y, const Real& z, const UInt& ic)
260  {
261  Real value (0.);
262  switch ( problemSolution )
263  {
264  case POISSON_POLYNOMIAL: // Poisson problem, polynomial solution
265  return -diffusionCoeff * 6.; // Polynomial solution
266  case POISSON_TRIGONOMETRIC: // Poisson problem, trigonometric solution
267  return 12 * Pi * Pi * std::sin (2 * Pi * x) * std::sin (2 * Pi * y) * std::sin (2 * Pi * z); // Trigonometric solution
268  case ADR_UNSTEADY_POLYNOMIAL: // ADR problem, polynomial solution
269  value += 2 * t;
270  case ADR_STEADY_POLYNOMIAL: // ADR problem, polynomial solution
271  value += -diffusionCoeff * 6 + reactionCoeff * u_ex (t, x, y, z, ic);
272  for (UInt icomp = 1; icomp <= nDimensions; ++icomp)
273  {
274  value += fAdvection (t, x, y, z, icomp) * grad_ex (icomp, t, x, y, z, ic);
275  }
276  return value;
277  default:
278  return 0;
279  }
280  }
281  // The advection term
282  static Real fAdvection (const Real& /*t*/, const Real& x, const Real& y, const Real& z, const UInt& icomp)
283  {
284  switch (icomp)
285  {
286  case 1:
287  switch ( problemSolution )
288  {
289  case POISSON_POLYNOMIAL: // Poisson problem, polynomial solution
290  case POISSON_TRIGONOMETRIC: // Poisson problem, trigonometric solution
291  return 0.;
292  case ADR_STEADY_POLYNOMIAL: // ADR problem, polynomial solution
293  case ADR_UNSTEADY_POLYNOMIAL: // ADR problem, polynomial solution
294  return 0.; //-0.1*x;
295  }
296  case 2:
297  switch ( problemSolution )
298  {
299  case POISSON_POLYNOMIAL: // Poisson problem, polynomial solution
300  case POISSON_TRIGONOMETRIC: // Poisson problem, trigonometric solution
301  return 0.;
302  case ADR_STEADY_POLYNOMIAL: // ADR problem, polynomial solution
303  case ADR_UNSTEADY_POLYNOMIAL: // ADR problem, polynomial solution
304  return 0.; //-0.1*y;
305  }
306  case 3:
307  switch ( problemSolution )
308  {
309  case POISSON_POLYNOMIAL: // Poisson problem, polynomial solution
310  case POISSON_TRIGONOMETRIC: // Poisson problem, trigonometric solution
311  return 0.;
312  case ADR_STEADY_POLYNOMIAL: // ADR problem, polynomial solution
313  case ADR_UNSTEADY_POLYNOMIAL: // ADR problem, polynomial solution
314  return 0.; //-0.1*z;
315  }
316  default:
317  return 0.;
318  }
319  }
320  // The normal derivative to the domain boundary
321  // BEWARE: this is case dependent! (works for cubic domains)
322  static Real fNeumann (const Real& t, const Real& x, const Real& y, const Real& z, const UInt& ic)
323  {
324  Real n[3] = {0, 0, 0};
325  Real out = 0;
326  if ( x == xMin )
327  {
328  n[0] = -1.;
329  }
330  else if ( x == xMax )
331  {
332  n[0] = 1.;
333  }
334  else if ( y == yMin )
335  {
336  n[1] = -1.;
337  }
338  else if ( y == yMax )
339  {
340  n[1] = 1.;
341  }
342  else if ( z == zMin )
343  {
344  n[2] = -1.;
345  }
346  else if ( z == zMax )
347  {
348  n[2] = 1.;
349  }
350  else
351  {
352  std::cout << "strange point: x=" << x << " y=" << y << " z=" << z
353  << std::endl;
354  }
355 
356  for (UInt icomp = 0; icomp < nDimensions; ++icomp) //mu grad_u n
357  {
358  out += diffusionCoeff * grad_ex (icomp + 1, t, x, y, z, ic) * n[icomp];
359  }
360  return out;
361  }
362  // Robin boundary condition compatible with the exact solution
363  static Real fRobin (const Real& t, const Real& x, const Real& y, const Real& z, const UInt& ic)
364  {
365  return fNeumann (t, x, y, z, ic) + coefRobin (t, x, y, z, ic) * u_ex (t, x, y, z, ic);
366  }
367  // A generic diffusion coefficient
368  //static Real diffusionCoeff(const Real& t, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const UInt& /*ic*/)
369  //{
370  // return t*t;
371  //}
372  // Initialization of the private members
373  static void setup ( const GetPot& dataFile, const std::string& section = "poisson" )
374  {
375  // We want a slash dividing the data file section from the variable name but only
376  // when not provided by the user or when not looking in the root of the data file
377  std::string corrected_section ( section );
378  if ( ( ! section.empty() ) && ( section[section.length() - 1] != '/' ) )
379  {
380  corrected_section = section + '/';
381  }
382 
383  xMin = dataFile ( (corrected_section + "problem/xMin").data(), -1. ) ;
384  xMax = dataFile ( (corrected_section + "problem/xMax").data(), 1. ) ;
385  yMin = dataFile ( (corrected_section + "problem/yMin").data(), -1. ) ;
386  yMax = dataFile ( (corrected_section + "problem/yMax").data(), 1. ) ;
387  zMin = dataFile ( (corrected_section + "problem/zMin").data(), -1. ) ;
388  zMax = dataFile ( (corrected_section + "problem/zMax").data(), 1. ) ;
389 
390  diffusionCoeff = dataFile ( (corrected_section + "physics/diffusionCoefficient").data(), 1. ) ;
391  reactionCoeff = dataFile ( (corrected_section + "physics/reactionCoefficient").data(), 1. ) ;
392  }
393  // Output on screen
394  static void showMe ( std::ostream& c = std::cout )
395  {
396  c << "\n*** Analytical solution: values for user-defined data\n";
397 
398  c << "\n[/problem]" << std::endl;
399  c << "xMin\t\t= " << xMin << std::endl;
400  c << "xMax\t\t= " << xMax << std::endl;
401  c << "yMin\t\t= " << yMin << std::endl;
402  c << "yMax\t\t= " << yMax << std::endl;
403  c << "zMin\t\t= " << zMin << std::endl;
404  c << "zMax\t\t= " << zMax << std::endl;
405  c << "solutionType\t\t= " << problemSolution << std::endl;
406 
407  c << "\n[/physics]" << std::endl;
408  c << "diffusionCoeff\t\t= " << diffusionCoeff << std::endl;
409  c << "reactionCoeff\t\t= " << reactionCoeff << std::endl;
410 
411  c << std::endl;
412  }
413 
414  // The dimensions of the domain
415  static Real xMin;
416  static Real xMax;
417  static Real yMin;
418  static Real yMax;
419  static Real zMin;
420  static Real zMax;
421  static Real diffusionCoeff;
422  static Real reactionCoeff;
423 
424  static ADRProblemSolution problemSolution;
425 
426 };
427 
428 // Declaration and default initialization of the static variables
429 Real AnalyticalSol::xMin (0.);
430 Real AnalyticalSol::xMax (1.);
431 Real AnalyticalSol::yMin (0.);
432 Real AnalyticalSol::yMax (1.);
433 Real AnalyticalSol::zMin (0.);
434 Real AnalyticalSol::zMax (1.);
435 Real AnalyticalSol::diffusionCoeff (1.);
436 Real AnalyticalSol::reactionCoeff (1.);
437 
438 ADRProblemSolution AnalyticalSol::problemSolution (POISSON_POLYNOMIAL);
439 
440 #endif
441 
442 
443 } //namespace LifeV
444 
445 #endif
LifeV::ADRProblemSolution
ADRProblemSolution
Definition: analyticalSol.hpp:38
LifeV::ADR_UNSTEADY_POLYNOMIAL
A polynomial solution to a general ADR problem.
Definition: analyticalSol.hpp:43
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::POISSON_TRIGONOMETRIC
A trigonometric solution to the Poisson problem.
Definition: analyticalSol.hpp:41
LifeV::ID
uint32_type ID
IDs.
Definition: LifeV.hpp:194
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::ADR_STEADY_POLYNOMIAL
A polynomial solution to a general ADR problem.
Definition: analyticalSol.hpp:42
LifeV::coefRobin
Real coefRobin(const Real &, const Real &, const Real &, const Real &, const ID &)
Definition: analyticalSol.hpp:48
LifeV::POISSON_POLYNOMIAL
A polynomial solution to the Poisson problem.
Definition: analyticalSol.hpp:40