35 #include "libmesh/quadrature.h"
36 #include "libmesh/boundary_info.h"
60 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
61 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
72 #ifdef GRINS_USE_GRVY_TIMERS
73 this->_timer->BeginTimer(
"AveragedTurbine::element_time_derivative");
77 const std::vector<libMesh::Real> &JxW =
78 context.get_element_fe(this->_flow_vars.u_var())->get_JxW();
81 const std::vector<std::vector<libMesh::Real> >& u_phi =
82 context.get_element_fe(this->_flow_vars.u_var())->get_phi();
84 const std::vector<libMesh::Point>& u_qpoint =
85 context.get_element_fe(this->_flow_vars.u_var())->get_xyz();
88 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u_var()).size();
91 libMesh::DenseSubMatrix<libMesh::Number> &Kuu = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.u_var());
92 libMesh::DenseSubMatrix<libMesh::Number> &Kuv = context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.v_var());
93 libMesh::DenseSubMatrix<libMesh::Number> &Kvu = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.u_var());
94 libMesh::DenseSubMatrix<libMesh::Number> &Kvv = context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.v_var());
96 libMesh::DenseSubMatrix<libMesh::Number> &Kus =
97 context.get_elem_jacobian(this->_flow_vars.u_var(),
98 this->fan_speed_var());
99 libMesh::DenseSubMatrix<libMesh::Number> &Ksu =
100 context.get_elem_jacobian(this->fan_speed_var(),
101 this->_flow_vars.u_var());
102 libMesh::DenseSubMatrix<libMesh::Number> &Kvs =
103 context.get_elem_jacobian(this->_flow_vars.v_var(),
104 this->fan_speed_var());
105 libMesh::DenseSubMatrix<libMesh::Number> &Ksv =
106 context.get_elem_jacobian(this->fan_speed_var(),
107 this->_flow_vars.v_var());
108 libMesh::DenseSubMatrix<libMesh::Number> &Kss =
109 context.get_elem_jacobian(this->fan_speed_var(),
110 this->fan_speed_var());
112 libMesh::DenseSubMatrix<libMesh::Number>* Kwu = NULL;
113 libMesh::DenseSubMatrix<libMesh::Number>* Kwv = NULL;
114 libMesh::DenseSubMatrix<libMesh::Number>* Kww = NULL;
115 libMesh::DenseSubMatrix<libMesh::Number>* Kuw = NULL;
116 libMesh::DenseSubMatrix<libMesh::Number>* Kvw = NULL;
118 libMesh::DenseSubMatrix<libMesh::Number>* Ksw = NULL;
119 libMesh::DenseSubMatrix<libMesh::Number>* Kws = NULL;
121 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u_var());
122 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v_var());
123 libMesh::DenseSubVector<libMesh::Number>* Fw = NULL;
125 libMesh::DenseSubVector<libMesh::Number> &Fs = context.get_elem_residual(this->fan_speed_var());
127 if( this->_dim == 3 )
129 Kuw = &context.get_elem_jacobian(this->_flow_vars.u_var(), this->_flow_vars.w_var());
130 Kvw = &context.get_elem_jacobian(this->_flow_vars.v_var(), this->_flow_vars.w_var());
132 Kwu = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.u_var());
133 Kwv = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.v_var());
134 Kww = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->_flow_vars.w_var());
135 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
137 Ksw = &context.get_elem_jacobian(this->fan_speed_var(), this->_flow_vars.w_var());
138 Kws = &context.get_elem_jacobian(this->_flow_vars.w_var(), this->fan_speed_var());
140 Fw = &context.get_elem_residual(this->_flow_vars.w_var());
143 unsigned int n_qpoints = context.get_element_qrule().n_points();
145 for (
unsigned int qp=0; qp != n_qpoints; qp++)
148 libMesh::Number u, v, s;
149 u = context.interior_value(this->_flow_vars.u_var(), qp);
150 v = context.interior_value(this->_flow_vars.v_var(), qp);
151 s = context.interior_value(this->fan_speed_var(), qp);
153 libMesh::NumberVectorValue U(u,v);
155 U(2) = context.interior_value(this->_flow_vars.w_var(), qp);
157 libMesh::NumberVectorValue U_B_1;
158 libMesh::NumberVectorValue F;
159 libMesh::NumberTensorValue dFdU;
160 libMesh::NumberTensorValue* dFdU_ptr =
161 compute_jacobian ? &dFdU : NULL;
162 libMesh::NumberVectorValue dFds;
163 libMesh::NumberVectorValue* dFds_ptr =
164 compute_jacobian ? &dFds : NULL;
165 if (!this->compute_force(u_qpoint[qp], context.time, U, s,
166 U_B_1, F, dFdU_ptr, dFds_ptr))
169 libMesh::Real jac = JxW[qp];
177 Fs(0) -= U_B_1 * F * jac;
179 if (compute_jacobian)
181 Kss(0,0) -= U_B_1 * dFds * jac;
183 for (
unsigned int j=0; j != n_u_dofs; j++)
185 libMesh::Real jac_j = JxW[qp] * u_phi[j][qp];
187 for (
unsigned int d=0; d != 3; ++d)
189 Ksu(0,j) -= jac_j * U_B_1(d) * dFdU(d,0);
190 Ksv(0,j) -= jac_j * U_B_1(d) * dFdU(d,1);
195 for (
unsigned int d=0; d != 3; ++d)
196 (*Ksw)(0,j) -= jac_j * U_B_1(d) * dFdU(d,2);
202 for (
unsigned int i=0; i != n_u_dofs; i++)
204 const libMesh::Number jac_i = jac * u_phi[i][qp];
209 if( this->_dim == 3 )
210 (*Fw)(i) += F(2)*jac_i;
212 if( compute_jacobian )
214 Kus(i,0) += dFds(0) * jac_i;
215 Kvs(i,0) += dFds(1) * jac_i;
216 if( this->_dim == 3 )
217 (*Kws)(i,0) += dFds(2) * jac_i;
219 for (
unsigned int j=0; j != n_u_dofs; j++)
221 const libMesh::Number jac_ij = jac_i * u_phi[j][qp];
222 Kuu(i,j) += jac_ij * dFdU(0,0);
223 Kuv(i,j) += jac_ij * dFdU(0,1);
224 Kvu(i,j) += jac_ij * dFdU(1,0);
225 Kvv(i,j) += jac_ij * dFdU(1,1);
227 if( this->_dim == 3 )
229 (*Kuw)(i,j) += jac_ij * dFdU(0,2);
230 (*Kvw)(i,j) += jac_ij * dFdU(1,2);
231 (*Kwu)(i,j) += jac_ij * dFdU(2,0);
232 (*Kwv)(i,j) += jac_ij * dFdU(2,1);
233 (*Kww)(i,j) += jac_ij * dFdU(2,2);
241 #ifdef GRINS_USE_GRVY_TIMERS
242 this->_timer->EndTimer(
"AveragedTurbine::element_time_derivative");
255 libMesh::DenseSubMatrix<libMesh::Number> &Kss =
256 context.get_elem_jacobian(this->fan_speed_var(), this->fan_speed_var());
258 libMesh::DenseSubVector<libMesh::Number> &Fs =
259 context.get_elem_residual(this->fan_speed_var());
261 const std::vector<libMesh::dof_id_type>& dof_indices =
262 context.get_dof_indices(this->fan_speed_var());
264 const libMesh::Number fan_speed =
265 context.get_system().current_solution(dof_indices[0]);
267 const libMesh::Number output_torque =
268 (*this->torque_function)(libMesh::Point(0), fan_speed);
270 Fs(0) += output_torque;
272 if (compute_jacobian)
275 const libMesh::Number epsilon = 1e-6;
276 const libMesh::Number output_torque_deriv =
277 ((*this->torque_function)(libMesh::Point(0), fan_speed+epsilon) -
278 (*this->torque_function)(libMesh::Point(0), fan_speed-epsilon)) / (2*epsilon);
280 Kss(0,0) += output_torque_deriv * context.get_elem_solution_derivative();
293 libMesh::DenseSubMatrix<libMesh::Number> &Kss =
294 context.get_elem_jacobian(this->fan_speed_var(), this->fan_speed_var());
296 libMesh::DenseSubVector<libMesh::Number> &Fs =
297 context.get_elem_residual(this->fan_speed_var());
299 const libMesh::DenseSubVector<libMesh::Number> &Us =
300 context.get_elem_solution_rate(this->_fan_speed_var);
302 const libMesh::Number& fan_speed = Us(0);
304 Fs(0) -= this->moment_of_inertia * fan_speed;
306 if (compute_jacobian)
308 Kss(0,0) -= this->moment_of_inertia * context.get_elem_solution_rate_derivative();
GRINS::ICHandlingBase * _ic_handler
virtual void nonlocal_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for scalar variables.
Base class for reading and handling initial conditions for physics classes.
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
virtual void nonlocal_mass_residual(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Mass matrix part(s) for scalar variables.
INSTANTIATE_INC_NS_SUBCLASS(AveragedTurbine)
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.