GRINS-0.8.0
velocity_penalty_base.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // GRINS - General Reacting Incompressible Navier-Stokes
5 //
6 // Copyright (C) 2014-2017 Paul T. Bauman, Roy H. Stogner
7 // Copyright (C) 2010-2013 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 
26 // This class
28 
29 // GRINS
31 
32 // libMesh
33 #include "libmesh/parsed_fem_function.h"
34 
35 namespace GRINS
36 {
37 
38  template<class Mu>
39  VelocityPenaltyBase<Mu>::VelocityPenaltyBase( const std::string& physics_name, const GetPot& input )
40  : IncompressibleNavierStokesBase<Mu>(physics_name,
41  PhysicsNaming::incompressible_navier_stokes(), /* "core" Physics name */
42  input),
43  base_physics_name("VelocityPenalty"),
44  _quadratic_scaling(false),
45  _input(input)
46  {
49  base_physics_name.push_back('2');
50 
53  base_physics_name.push_back('3');
54 
55  this->read_input_options(input);
56  }
57 
58  template<class Mu>
59  void VelocityPenaltyBase<Mu>::read_input_options( const GetPot& input )
60  {
61  _quadratic_scaling =
62  input("Physics/"+base_physics_name+"/quadratic_scaling", false);
63  }
64 
65  template<class Mu>
66  void VelocityPenaltyBase<Mu>::set_time_evolving_vars ( libMesh::FEMSystem* system)
67  {
70  this->normal_vector_function.reset(nvf);
71 
72  this->set_parameter(*nvf, _input,
73  "Physics/"+base_physics_name+"/penalty_function",
74  this->zero_vector_function);
75 
78  this->base_velocity_function.reset(bvf);
79 
80  this->set_parameter(*bvf, _input,
81  "Physics/"+base_physics_name+"/base_velocity",
82  this->zero_vector_function);
83 
84  if ((nvf->expression() == this->zero_vector_function) &&
85  (bvf->expression() == this->zero_vector_function))
86  std::cout << "Warning! Zero VelocityPenalty specified!" << std::endl;
87  }
88 
89 
90  template<class Mu>
92  ( const libMesh::Point& point,
93  const AssemblyContext& context,
94  const libMesh::NumberVectorValue& U,
95  libMesh::NumberVectorValue& F,
96  libMesh::NumberTensorValue *dFdU)
97  {
98  // Velocity discrepancy (current velocity minus base velocity)
99  // normal to constraint plane, scaled by constraint penalty
100  // value
101  libMesh::DenseVector<libMesh::Number> output_vec(3);
102 
103  (*this->normal_vector_function)(context, point, context.time, output_vec);
104 
105  libMesh::NumberVectorValue U_N(output_vec(0),
106  output_vec(1),
107  output_vec(2));
108 
109  (*this->base_velocity_function)(context, point, context.time, output_vec);
110 
111  const libMesh::NumberVectorValue U_B(output_vec(0),
112  output_vec(1),
113  output_vec(2));
114 
115  const libMesh::NumberVectorValue U_Rel = U-U_B;
116 
117  // Old code
118  // const libMesh::NumberVectorValue F1 = (U_Rel*U_N)*U_N; //
119 
120  // With correct sign and more natural normalization
121  const libMesh::Number U_N_mag = std::sqrt(U_N*U_N);
122 
123  if (!U_N_mag)
124  return false;
125 
126  const libMesh::NumberVectorValue U_N_unit = U_N/U_N_mag;
127 
128  F = -(U_Rel*U_N)*U_N_unit;
129 
130  if (dFdU)
131  for (unsigned int i=0; i != 3; ++i)
132  for (unsigned int j=0; j != 3; ++j)
133  (*dFdU)(i,j) = -(U_N(j))*U_N_unit(i);
134 
135  // With quadratic scaling
136  if (_quadratic_scaling)
137  {
138  const libMesh::Number U_Rel_mag = std::sqrt(U_Rel * U_Rel);
139 
140  // Modify dFdU first so as to reuse the old value of F
141  if (dFdU)
142  {
143  // dU_Rel/dU = I
144  // d(U_Rel*U_Rel)/dU = 2*U_Rel
145  // d|U_Rel|/dU = U_Rel/|U_Rel|
146 
147  const libMesh::NumberVectorValue U_Rel_unit = U_Rel/U_Rel_mag;
148 
149  (*dFdU) *= U_Rel_mag;
150 
151  if (U_Rel_mag)
152  for (unsigned int i=0; i != 3; ++i)
153  for (unsigned int j=0; j != 3; ++j)
154  (*dFdU)(i,j) += F(i)*U_Rel_unit(j);
155  }
156 
157  F *= U_Rel_mag;
158  }
159 
160  // With correction term to avoid doing work on flow
161  bool do_correction = false;
162  if (do_correction)
163  {
164  if (dFdU)
165  libmesh_not_implemented();
166 
167  const libMesh::Number U_Rel_mag_sq = U_Rel * U_Rel;
168  if (U_Rel_mag_sq)
169  {
170  F -= (U_Rel*F)*U_Rel/U_Rel_mag_sq;
171  }
172  }
173  return true;
174  }
175 
176 } // namespace GRINS
177 
178 // Instantiate
179 INSTANTIATE_INC_NS_SUBCLASS(VelocityPenaltyBase);
Physics class for Incompressible Navier-Stokes.
static PhysicsName velocity_penalty3_adjoint_stab()
bool compute_force(const libMesh::Point &point, const AssemblyContext &context, const libMesh::NumberVectorValue &U, libMesh::NumberVectorValue &F, libMesh::NumberTensorValue *dFdU=NULL)
GRINS namespace.
static PhysicsName velocity_penalty2()
const PhysicsName _physics_name
Name of the physics object. Used for reading physics specific inputs.
Definition: physics.h:256
INSTANTIATE_INC_NS_SUBCLASS(VelocityPenaltyBase)
void read_input_options(const GetPot &input)
Read options from GetPot input file.
void set_time_evolving_vars(libMesh::FEMSystem *system)
Sets velocity variables to be time-evolving.
static PhysicsName velocity_penalty2_adjoint_stab()
static PhysicsName velocity_penalty3()

Generated on Tue Dec 19 2017 12:47:28 for GRINS-0.8.0 by  doxygen 1.8.9.1