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