LifeV
RossEthierSteinmanInc.cpp
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 /*!
28  @file Ross-Ethier Steinmann Analytical Solution
29  @brief A short description of the test content
30 
31  @author Christoph Winkelmann <christoph.winkelmann@epfl.ch>
32  @contributor Mauro Perego <mauro@mathcs.emory.edu>
33  @contributor Umberto Villa <uvilla@emory.edu>
34  @maintainer Umberto Villa <uvilla@emory.edu>
35 
36  @date 2004-11-12
37 
38  Implementation.
39 */
40 
41 #include <lifev/core/LifeV.hpp>
42 
43 #include <lifev/navier_stokes/function/RossEthierSteinmanInc.hpp>
44 
45 namespace LifeV
46 {
47 
48 
49 Real RossEthierSteinmanUnsteadyInc::uexact ( const Real& t,
50  const Real& x,
51  const Real& y,
52  const Real& z,
53  const ID& i)
54 {
55  Real e = std::exp (2.* (S_a * S_a + S_d * S_d + S_a * S_d) * S_nu * t);
56  switch (i)
57  {
58  case 0:
59  return e * ( S_d * std::exp (S_a * (x - z) + S_d * (y - z) ) - S_a * std::exp (S_a * (z - y) + S_d * (x - y) ) );
60  case 1:
61  return e * ( S_d * std::exp (S_a * (y - x) + S_d * (z - x) ) - S_a * std::exp (S_a * (x - z) + S_d * (y - z) ) );
62  case 2:
63  return e * ( S_d * std::exp (S_a * (z - y) + S_d * (x - y) ) - S_a * std::exp (S_a * (y - x) + S_d * (z - x) ) );
64  default:
65  exit (1);
66  }
67  return 1.;
68 }
69 
70 Real RossEthierSteinmanUnsteadyInc::uderexact ( const Real& t,
71  const Real& x,
72  const Real& y,
73  const Real& z,
74  const ID& i)
75 {
76 
77  if (i < 3)
78  {
79  return 2.* (S_a * S_a + S_d * S_d + S_a * S_d) * S_nu * xexact (t, x, y, z, i);
80  }
81  else
82  {
83  return 0.;
84  }
85 }
86 
87 Real RossEthierSteinmanUnsteadyInc::pexact ( const Real& t,
88  const Real& x,
89  const Real& y,
90  const Real& z,
91  const ID& /* i */ )
92 {
93  return S_rho * (S_a * S_a + S_d * S_d + S_a * S_d) * std::exp (4.* (S_a * S_a + S_d * S_d + S_a * S_d) * S_nu * t) *
94  (std::exp (S_a * (x - y) + S_d * (x - z) ) + std::exp (S_a * (y - z) + S_d * (y - x) ) + std::exp (S_a * (z - x) + S_d * (z - y) ) );
95 }
96 
97 Real RossEthierSteinmanUnsteadyInc::grad_u ( const UInt& icoor, const Real& t, const Real& x, const Real& y, const Real& z, const ID& i )
98 {
99  Real e = std::exp (2.* (S_a * S_a + S_d * S_d + S_a * S_d) * S_nu * t);
100  switch (icoor)
101  {
102  case 0: // u_x
103  switch (i)
104  {
105  case 0:
106  return S_a * S_d * e * ( std::exp (S_a * (x - z) + S_d * (y - z) ) - std::exp (S_a * (z - y) + S_d * (x - y) ) );
107  case 1:
108  return e * ( - (S_a + S_d) * S_d * std::exp (S_a * (y - x) + S_d * (z - x) ) - S_a * S_a * std::exp (S_a * (x - z) + S_d * (y - z) ) );
109  case 2:
110  return e * ( S_d * S_d * std::exp (S_a * (z - y) + S_d * (x - y) ) + S_a * (S_a + S_d) * std::exp (S_a * (y - x) + S_d * (z - x) ) );
111  default:
112  exit (1);
113  }
114  case 1: // u_y
115  switch (i)
116  {
117  case 0:
118  return e * ( S_d * S_d * std::exp (S_a * (x - z) + S_d * (y - z) ) + S_a * (S_a + S_d) * std::exp (S_a * (z - y) + S_d * (x - y) ) );
119  case 1:
120  return e * ( S_a * S_d * std::exp (S_a * (y - x) + S_d * (z - x) ) - S_a * S_d * std::exp (S_a * (x - z) + S_d * (y - z) ) );
121  case 2:
122  return e * ( -S_d * (S_a + S_d) * std::exp (S_a * (z - y) + S_d * (x - y) ) - S_a * S_a * std::exp (S_a * (y - x) + S_d * (z - x) ) );
123  default:
124  exit (1);
125  }
126  case 2:
127  switch (i)
128  {
129  case 0:
130  return e * ( -S_d * (S_a + S_d) * std::exp (S_a * (x - z) + S_d * (y - z) ) - S_a * S_a * std::exp (S_a * (z - y) + S_d * (x - y) ) );
131  case 1:
132  return e * ( S_d * S_d * std::exp (S_a * (y - x) + S_d * (z - x) ) + S_a * (S_a + S_d) * std::exp (S_a * (x - z) + S_d * (y - z) ) );
133  case 2:
134  return e * ( S_a * S_d * std::exp (S_a * (z - y) + S_d * (x - y) ) - S_a * S_d * std::exp (S_a * (y - x) + S_d * (z - x) ) );
135  default:
136  exit (1);
137  }
138  default:
139  exit (1);
140  }
141  return 1.;
142 }
143 
144 Real RossEthierSteinmanUnsteadyInc::f ( const Real& /* t */,
145  const Real& /* x */,
146  const Real& /* y */,
147  const Real& /* z */,
148  const ID& /* i */ )
149 {
150  return 0.;
151 }
152 
153 Real RossEthierSteinmanUnsteadyInc::xexact ( const Real& t,
154  const Real& x,
155  const Real& y,
156  const Real& z,
157  const ID& i )
158 {
159  switch (i)
160  {
161  case 0:
162  case 1:
163  case 2:
164  return uexact (t, x, y, z, i);
165  case 3:
166  return pexact (t, x, y, z, 1);
167  default:
168  exit (1);
169  }
170  return 1.;
171 }
172 
173 
174 
175 // Initial velocity
176 Real RossEthierSteinmanUnsteadyInc::x0 ( const Real& t, const Real& x, const Real& y,
177  const Real& z, const ID& i )
178 {
179  return xexact (t, x, y, z, i);
180 }
181 
182 
183 
184 //we suppose that the problem geometry is the cube [-1,1]x[-1,1]x[-1,1].
185 Real RossEthierSteinmanUnsteadyInc::fNeumann ( const Real& t,
186  const Real& x,
187  const Real& y,
188  const Real& z,
189  const ID& i )
190 {
191  Real n[3] = {0., 0., 0.};
192  Real out = 0.;
193  if ( x == -1. )
194  {
195  n[0] = -1.;
196  }
197  else if ( x == 1. )
198  {
199  n[0] = 1.;
200  }
201  else if ( y == -1. )
202  {
203  n[1] = -1.;
204  }
205  else if ( y == 1. )
206  {
207  n[1] = 1.;
208  }
209  else if ( z == -1. )
210  {
211  n[2] = -1.;
212  }
213  else if ( z == 1. )
214  {
215  n[2] = 1.;
216  }
217  else
218  {
219  std::cout << "strange point: x=" << x << " y=" << y << " z=" << z
220  << std::endl;
221  }
222 
223  for (UInt k = 0; k < 3; k++) //mu grad_u n
224  {
225  out += S_mu * grad_u (k, t, x, y, z, i) * n[k];
226  }
227 
228  if (S_flagStrain)
229  for (UInt k = 0; k < 3; k++) //mu grad_u^T n
230  {
231  out += S_mu * grad_u (i, t, x, y, z, k) * n[k];
232  }
233 
234  out -= pexact (t, x, y, z, i) * n[i]; //grad_p n
235  return out;
236 }
237 
238 
239 void RossEthierSteinmanUnsteadyInc::setParamsFromGetPot ( const GetPot& dataFile )
240 {
241  S_a = dataFile ( "fluid/problem/a", 1. ) ;
242  S_d = dataFile ( "fluid/problem/d", 1. ) ;
243  S_mu = dataFile ( "fluid/physics/viscosity", 1. );
244  S_rho = dataFile ( "fluid/physics/density", 1. );
245  S_nu = S_mu / S_rho;
246  S_flagStrain = dataFile ( "fluid/physics/flag_strain", 0 );
247 }
248 void RossEthierSteinmanUnsteadyInc::setA (const Real& aValue)
249 {
250  S_a = aValue;
251 }
252 void RossEthierSteinmanUnsteadyInc::setD (const Real& dValue)
253 {
254  S_d = dValue;
255 }
256 void RossEthierSteinmanUnsteadyInc::setViscosity (const Real& mu)
257 {
258  S_mu = mu;
259  S_nu = S_mu / S_rho;
260 }
261 void RossEthierSteinmanUnsteadyInc::setDensity (const Real& rho)
262 {
263  S_rho = rho;
264  S_nu = S_mu / S_rho;
265 }
266 void RossEthierSteinmanUnsteadyInc::setFlagStrain (const Int& flagValue)
267 {
268  S_flagStrain = flagValue;
269 }
270 
271 Real RossEthierSteinmanUnsteadyInc::S_a;
272 Real RossEthierSteinmanUnsteadyInc::S_d;
273 Real RossEthierSteinmanUnsteadyInc::S_mu;
274 Real RossEthierSteinmanUnsteadyInc::S_rho;
275 Real RossEthierSteinmanUnsteadyInc::S_nu;
276 Int RossEthierSteinmanUnsteadyInc::S_flagStrain;
277 
278 } // namespace LifeV
GetPot::operator()
double operator()(const char *VarName, const double &Default) const
Definition: GetPot.hpp:2033
LifeV::Int
int32_type Int
Generic integer data.
Definition: LifeV.hpp:188
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::ID
uint32_type ID
IDs.
Definition: LifeV.hpp:194
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
GetPot::operator()
int operator()(const char *VarName, int Default) const
Definition: GetPot.hpp:2021
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191