www.mooseframework.org
RichardsSUPGstandard.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "RichardsSUPGstandard.h"
9 #include "libmesh/utility.h"
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<RichardsSUPG>();
16  params.addRequiredRangeCheckedParam<Real>(
17  "p_SUPG",
18  "p_SUPG > 0",
19  "SUPG pressure. This parameter controls the strength of the "
20  "upwinding. This parameter must be positive. If you need to track "
21  "advancing fronts in a problem, then set to less than your expected "
22  "range of pressures in your unsaturated zone. Otherwise, set "
23  "larger, and then minimal upwinding will occur and convergence will "
24  "typically be good. If you need no SUPG, it is more efficient not "
25  "to use this UserObject.");
26  params.addClassDescription(
27  "Standard SUPG relationships for Richards flow based on Appendix A of TJR Hughes, M "
28  "Mallet and A Mizukami ``A new finite element formulation for computational fluid dynamics:: "
29  "II. Behond SUPG'' Computer Methods in Applied Mechanics and Engineering 54 (1986) 341--355");
30  return params;
31 }
32 
33 RichardsSUPGstandard::RichardsSUPGstandard(const InputParameters & parameters)
34  : RichardsSUPG(parameters), _p_SUPG(getParam<Real>("p_SUPG"))
35 {
36 }
37 
38 RealVectorValue
39 RichardsSUPGstandard::velSUPG(RealTensorValue perm,
40  RealVectorValue gradp,
41  Real density,
42  RealVectorValue gravity) const
43 {
44  return -perm * (gradp - density * gravity); // points in direction of info propagation
45 }
46 
47 RealTensorValue
48 RichardsSUPGstandard::dvelSUPG_dgradp(RealTensorValue perm) const
49 {
50  return -perm;
51 }
52 
53 RealVectorValue
54 RichardsSUPGstandard::dvelSUPG_dp(RealTensorValue perm,
55  Real density_prime,
56  RealVectorValue gravity) const
57 {
58  return perm * density_prime * gravity;
59 }
60 
61 Real
63 {
64  if (alpha >= 5.0 || alpha <= -5.0)
65  return ((alpha > 0.0) ? 1.0 : -1.0) - 1.0 / alpha; // prevents overflows
66  else if (alpha == 0)
67  return 0.0;
68  return 1.0 / std::tanh(alpha) - 1.0 / alpha;
69 }
70 
71 Real
73 {
74  if (alpha >= 5.0 || alpha <= -5.0)
75  return 1.0 / (alpha * alpha); // prevents overflows
76  else if (alpha == 0)
77  return 1.0 / 3.0;
78  return 1.0 - Utility::pow<2>(std::cosh(alpha) / std::sinh(alpha)) + 1.0 / (alpha * alpha);
79 }
80 
81 // the following is physically 2*velocity/element_length
82 RealVectorValue
83 RichardsSUPGstandard::bb(RealVectorValue vel,
84  int dimen,
85  RealVectorValue xi_prime,
86  RealVectorValue eta_prime,
87  RealVectorValue zeta_prime) const
88 {
89  RealVectorValue b;
90  b(0) = xi_prime * vel;
91  if (dimen >= 2)
92  b(1) = eta_prime * vel;
93  if (dimen == 3)
94  b(2) = zeta_prime * vel;
95  return b;
96 }
97 
98 // following is d(bb*bb)/d(gradp)
99 RealVectorValue
101  RealTensorValue dvel_dgradp,
102  RealVectorValue xi_prime,
103  RealVectorValue eta_prime,
104  RealVectorValue zeta_prime) const
105 {
106  // if dvel_dgradp is symmetric so transpose is probably a waste of time
107  return 2.0 * ((xi_prime * vel) * (dvel_dgradp.transpose() * xi_prime) +
108  (eta_prime * vel) * (dvel_dgradp.transpose() * eta_prime) +
109  (zeta_prime * vel) * (dvel_dgradp.transpose() * zeta_prime));
110 }
111 
112 // following is d(bb*bb)/d(p)
113 Real
114 RichardsSUPGstandard::dbb2_dp(RealVectorValue vel,
115  RealVectorValue dvel_dp,
116  RealVectorValue xi_prime,
117  RealVectorValue eta_prime,
118  RealVectorValue zeta_prime) const
119 {
120  return 2.0 *
121  ((xi_prime * vel) * (dvel_dp * xi_prime) + (eta_prime * vel) * (dvel_dp * eta_prime) +
122  (zeta_prime * vel) * (dvel_dp * zeta_prime));
123 }
124 
125 Real
126 RichardsSUPGstandard::tauSUPG(RealVectorValue vel, Real traceperm, RealVectorValue b) const
127 {
128  // vel = velocity, b = bb
129  Real norm_v = vel.norm();
130  Real norm_b = b.norm(); // Hughes et al investigate infinity-norm and 2-norm. i just use 2-norm
131  // here. norm_b ~ 2|a|/ele_length_in_direction_of_a
132 
133  if (norm_b == 0)
134  return 0.0; // Only way for norm_b=0 is for zero ele size, or vel=0. Either way we don't have
135  // to upwind.
136 
137  Real h = 2.0 * norm_v / norm_b; // h is a measure of the element length in the "a" direction
138  Real alpha = 0.5 * norm_v * h / (traceperm * _p_SUPG); // this is the Peclet number
139 
140  const Real xi_tilde = RichardsSUPGstandard::cosh_relation(alpha);
141 
142  return xi_tilde / norm_b;
143 }
144 
145 RealVectorValue
147  RealTensorValue dvel_dgradp,
148  Real traceperm,
149  RealVectorValue b,
150  RealVectorValue db2_dgradp) const
151 {
152  Real norm_vel = vel.norm();
153  if (norm_vel == 0)
154  return RealVectorValue();
155 
156  RealVectorValue norm_vel_dgradp(dvel_dgradp * vel / norm_vel);
157 
158  Real norm_b = b.norm();
159  if (norm_b == 0)
160  return RealVectorValue();
161  RealVectorValue norm_b_dgradp = db2_dgradp / 2.0 / norm_b;
162 
163  Real h = 2 * norm_vel / norm_b; // h is a measure of the element length in the "a" direction
164  RealVectorValue h_dgradp(2 * norm_vel_dgradp / norm_b -
165  2.0 * norm_vel * norm_b_dgradp / norm_b / norm_b);
166 
167  Real alpha = 0.5 * norm_vel * h / traceperm / _p_SUPG; // this is the Peclet number
168  RealVectorValue alpha_dgradp =
169  0.5 * (norm_vel_dgradp * h + norm_vel * h_dgradp) / traceperm / _p_SUPG;
170 
171  Real xi_tilde = RichardsSUPGstandard::cosh_relation(alpha);
172  Real xi_tilde_prime = RichardsSUPGstandard::cosh_relation_prime(alpha);
173  RealVectorValue xi_tilde_dgradp = xi_tilde_prime * alpha_dgradp;
174 
175  RealVectorValue tau_dgradp =
176  xi_tilde_dgradp / norm_b - xi_tilde * norm_b_dgradp / (norm_b * norm_b);
177 
178  return tau_dgradp;
179 }
180 
181 Real
183  RealVectorValue dvel_dp,
184  Real traceperm,
185  RealVectorValue b,
186  Real db2_dp) const
187 {
188  Real norm_vel = vel.norm();
189  if (norm_vel == 0.0)
190  return 0.0; // this deriv is not necessarily correct, but i can't see a better thing to do
191 
192  Real norm_vel_dp = dvel_dp * vel / norm_vel;
193 
194  Real norm_b = b.norm();
195  if (norm_b == 0)
196  return 0.0; // this deriv is not necessarily correct, but i can't see a better thing to do
197  Real norm_b_dp = db2_dp / (2.0 * norm_b);
198 
199  Real h = 2.0 * norm_vel / norm_b; // h is a measure of the element length in the "a" direction
200  Real h_dp = 2.0 * norm_vel_dp / norm_b - 2.0 * norm_vel * norm_b_dp / (norm_b * norm_b);
201 
202  Real alpha = 0.5 * norm_vel * h / (traceperm * _p_SUPG); // this is the Peclet number
203  Real alpha_dp = 0.5 * (norm_vel_dp * h + norm_vel * h_dp) / (traceperm * _p_SUPG);
204 
205  Real xi_tilde = RichardsSUPGstandard::cosh_relation(alpha);
206  Real xi_tilde_prime = RichardsSUPGstandard::cosh_relation_prime(alpha);
207  Real xi_tilde_dp = xi_tilde_prime * alpha_dp;
208 
209  // Real tau = xi_tilde/norm_b;
210  const Real tau_dp = xi_tilde_dp / norm_b - xi_tilde * norm_b_dp / (norm_b * norm_b);
211 
212  return tau_dp;
213 }
214 
215 bool
217 {
218  return false;
219 }
RichardsSUPGstandard(const InputParameters &parameters)
Real _p_SUPG
the SUPG pressure parameter This dictates how strong the SUPG is _p_SUPG large means only a little SU...
RealVectorValue dbb2_dgradp(RealVectorValue vel, RealTensorValue dvel_dgradp, RealVectorValue xi_prime, RealVectorValue eta_prime, RealVectorValue zeta_prime) const
derivative of bb*bb wrt gradient of porepressure
const std::string density
Definition: NS.h:15
InputParameters validParams< RichardsSUPGstandard >()
Real dtauSUPG_dp(RealVectorValue vel, RealVectorValue dvel_dp, Real traceperm, RealVectorValue b, Real db2_dp) const
derivative of tau wrt porepressure (keeping gradp fixed)
Real cosh_relation_prime(Real alpha) const
derivative of cosh_relation wrt alpha
base class for SUPG of the Richards equation You must override all the functions below with your spec...
Definition: RichardsSUPG.h:24
RealTensorValue dvelSUPG_dgradp(RealTensorValue perm) const
derivative of SUPG velocity wrt gradient of porepressure
RealVectorValue bb(RealVectorValue vel, int dimen, RealVectorValue xi_prime, RealVectorValue eta_prime, RealVectorValue zeta_prime) const
|bb| ~ 2*velocity/element_length
Real tauSUPG(RealVectorValue vel, Real traceperm, RealVectorValue b) const
The SUPG tau parameter.
bool SUPG_trivial() const
returns false in this case since everything is trivial
InputParameters validParams< RichardsSUPG >()
Definition: RichardsSUPG.C:14
RealVectorValue dtauSUPG_dgradp(RealVectorValue vel, RealTensorValue dvel_dgradp, Real traceperm, RealVectorValue b, RealVectorValue db2_dgradp) const
derivative of tau wrt gradient of porepressure
RealVectorValue velSUPG(RealTensorValue perm, RealVectorValue gradp, Real density, RealVectorValue gravity) const
SUPG velocity = -perm*(gradp - density*gravity) This points in direction of information propagation...
Real cosh_relation(Real alpha) const
cosh(alpha)/sinh(alpha) - 1/alpha, modified at extreme values of alpha to prevent overflows ...
RealVectorValue dvelSUPG_dp(RealTensorValue perm, Real density_prime, RealVectorValue gravity) const
derivative of SUPG velocity wrt porepressure (keeping gradp fixed)
Real dbb2_dp(RealVectorValue vel, RealVectorValue dvel_dp, RealVectorValue xi_prime, RealVectorValue eta_prime, RealVectorValue zeta_prime) const
derivative of bb*bb wrt porepressure