refineable_linearised_navier_stokes_elements.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Non-inline functions for the refineable linearised axisymmetric
27// Navier-Stokes elements
28
29// oomph-lib includes
31
32namespace oomph
33{
34 //=======================================================================
35 /// Compute the residuals for the refineable linearised axisymmetric
36 /// Navier--Stokes equations; flag=1(or 0): do (or don't) compute the
37 /// Jacobian as well.
38 //=======================================================================
41 Vector<double>& residuals,
42 DenseMatrix<double>& jacobian,
43 DenseMatrix<double>& mass_matrix,
44 unsigned flag)
45 {
46 // Get the time from the first node in the element
47 const double time = this->node_pt(0)->time_stepper_pt()->time();
48
49 // Determine number of nodes in the element
50 const unsigned n_node = nnode();
51
52 // Determine how many pressure values there are associated with
53 // a single pressure component
54 const unsigned n_pres = npres_linearised_nst();
55
56 const unsigned n_veloc = 4 * DIM;
57
58 // Get the nodal indices at which the velocity is stored
59 unsigned u_nodal_index[n_veloc];
60 for (unsigned i = 0; i < n_veloc; ++i)
61 {
62 u_nodal_index[i] = u_index_linearised_nst(i);
63 }
64
65 // Which nodal values represent the two pressure components?
66 // (Negative if pressure is not based on nodal interpolation).
67 Vector<int> p_index(2);
68 for (unsigned i = 0; i < 2; i++)
69 {
70 p_index[i] = this->p_index_linearised_nst(i);
71 }
72
73 // Local array of booleans that are true if the l-th pressure value is
74 // hanging (avoid repeated virtual function calls)
75 bool pressure_dof_is_hanging[n_pres];
76
77 // If the pressure is stored at a node
78 if (p_index[0] >= 0)
79 {
80 // Read out whether the pressure is hanging
81 for (unsigned l = 0; l < n_pres; ++l)
82 {
83 pressure_dof_is_hanging[l] =
84 pressure_node_pt(l)->is_hanging(p_index[0]);
85 }
86 }
87 // Otherwise the pressure is not stored at a node and so cannot hang
88 else
89 {
90 for (unsigned l = 0; l < n_pres; ++l)
91 {
92 pressure_dof_is_hanging[l] = false;
93 }
94 }
95
96
97 // Set up memory for the fluid shape and test functions
98 Shape psif(n_node), testf(n_node);
99 DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
100
101 // Set up memory for the pressure shape and test functions
102 Shape psip(n_pres), testp(n_pres);
103
104 // Determine number of integration points
105 const unsigned n_intpt = integral_pt()->nweight();
106
107 // Set up memory for the vector to hold local coordinates (two dimensions)
108 Vector<double> s(DIM);
109
110 // Get physical variables from the element
111 // (Reynolds number must be multiplied by the density ratio)
112 const double scaled_re = re() * density_ratio();
113 const double scaled_re_st = re_st() * density_ratio();
114 const double visc_ratio = viscosity_ratio();
115
116 const double eval_real = lambda();
117 const double eval_imag = omega();
118
119 const std::complex<double> eigenvalue(eval_real, eval_imag);
120
121 // Integers used to store the local equation numbers
122 int local_eqn = 0;
123
124 // Local storage for pointers to hang info objects
125 HangInfo* hang_info_pt = 0;
126
127 // Loop over the integration points
128 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
129 {
130 // Assign values of the local coordinates s
131 for (unsigned i = 0; i < DIM; i++)
132 {
133 s[i] = integral_pt()->knot(ipt, i);
134 }
135
136 // Get the integral weight
137 const double w = integral_pt()->weight(ipt);
138
139 // Calculate the fluid shape and test functions, and their derivatives
140 // w.r.t. the global coordinates
142 ipt, psif, dpsifdx, testf, dtestfdx);
143
144 // Calculate the pressure shape and test functions
145 pshape_linearised_nst(s, psip, testp);
146
147 // Premultiply the weights and the Jacobian of the mapping between
148 // local and global coordinates
149 const double W = w * J;
150
151 // Allocate storage for the position and the derivative of the
152 // mesh positions w.r.t. time
154 // Vector<double> mesh_velocity(2,0.0);
155
156 // Allocate storage for the velocity components (six of these)
157 // and their derivatives w.r.t. time
158 Vector<std::complex<double>> interpolated_u(DIM);
159 // Vector<double> dudt(6,0.0);
160 // Allocate storage for the eigen function normalisation
161 Vector<std::complex<double>> interpolated_u_normalisation(DIM);
162 for (unsigned i = 0; i < DIM; ++i)
163 {
164 interpolated_u[i].real(0.0);
165 interpolated_u[i].imag(0.0);
166 interpolated_u_normalisation[i].real(0.0);
167 interpolated_u_normalisation[i].imag(0.0);
168 }
169
170 // Allocate storage for the pressure components (two of these
171 std::complex<double> interpolated_p(0.0, 0.0);
172 std::complex<double> interpolated_p_normalisation(0.0, 0.0);
173
174 // Allocate storage for the derivatives of the velocity components
175 // w.r.t. global coordinates (r and z)
176 Vector<Vector<std::complex<double>>> interpolated_dudx(DIM);
177 for (unsigned i = 0; i < DIM; ++i)
178 {
179 interpolated_dudx[i].resize(DIM);
180 for (unsigned j = 0; j < DIM; ++j)
181 {
182 interpolated_dudx[i][j].real(0.0);
183 interpolated_dudx[i][j].imag(0.0);
184 }
185 }
186
187 // Calculate pressure at the integration point
188 // -------------------------------------------
189
190 // Loop over pressure degrees of freedom (associated with a single
191 // pressure component) in the element
192 for (unsigned l = 0; l < n_pres; l++)
193 {
194 // Cache the shape function
195 const double psip_ = psip(l);
196
197 // Get the complex nodal pressure values
198 const std::complex<double> p_value(this->p_linearised_nst(l, 0),
199 this->p_linearised_nst(l, 1));
200
201 // Add contribution
202 interpolated_p += p_value * psip_;
203
204 // Repeat for normalisation
205 const std::complex<double> p_norm_value(this->p_linearised_nst(l, 2),
206 this->p_linearised_nst(l, 3));
207 interpolated_p_normalisation += p_norm_value * psip_;
208 }
209 // End of loop over the pressure degrees of freedom in the element
210
211 // Calculate velocities and their derivatives at the integration point
212 // -------------------------------------------------------------------
213
214 // Loop over the element's nodes
215 for (unsigned l = 0; l < n_node; l++)
216 {
217 // Cache the shape function
218 const double psif_ = psif(l);
219
220 // Loop over the DIM coordinate directions
221 for (unsigned i = 0; i < DIM; i++)
222 {
223 interpolated_x[i] += this->nodal_position(l, i) * psif_;
224 }
225
226 // Loop over the DIM complex velocity components
227 for (unsigned i = 0; i < DIM; i++)
228 {
229 // Get the value
230 const std::complex<double> u_value(
231 this->nodal_value(l, u_nodal_index[2 * i + 0]),
232 this->nodal_value(l, u_nodal_index[2 * i + 1]));
233
234 // Add contribution
235 interpolated_u[i] += u_value * psif_;
236
237 // Add contribution to dudt
238 // dudt[i] += du_dt_linearised_nst(l,i)*psif_;
239
240 // Loop over two coordinate directions (for derivatives)
241 for (unsigned j = 0; j < DIM; j++)
242 {
243 interpolated_dudx[i][j] += u_value * dpsifdx(l, j);
244 }
245
246 // Interpolate the normalisation function
247 const std::complex<double> normalisation_value(
248 this->nodal_value(l, u_nodal_index[2 * (DIM + i)]),
249 this->nodal_value(l, u_nodal_index[2 * (DIM + i) + 1]));
250 interpolated_u_normalisation[i] += normalisation_value * psif_;
251 }
252 } // End of loop over the element's nodes
253
254 // Get the mesh velocity if ALE is enabled
255 /*if(!ALE_is_disabled)
256 {
257 // Loop over the element's nodes
258 for(unsigned l=0;l<n_node;l++)
259 {
260 // Loop over the two coordinate directions
261 for(unsigned i=0;i<2;i++)
262 {
263 mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif(l);
264 }
265 }
266 }*/
267
268 // Get velocities and their derivatives from base flow problem
269 // -----------------------------------------------------------
270
271 // Allocate storage for the velocity components of the base state
272 // solution (initialise to zero)
273 Vector<double> base_flow_u(DIM, 0.0);
274
275 // Get the user-defined base state solution velocity components
276 get_base_flow_u(time, ipt, interpolated_x, base_flow_u);
277
278 // Allocate storage for the derivatives of the base state solution's
279 // velocity components w.r.t. global coordinate (r and z)
280 // N.B. the derivatives of the base flow components w.r.t. the
281 // azimuthal coordinate direction (theta) are always zero since the
282 // base flow is axisymmetric
283 DenseMatrix<double> base_flow_dudx(DIM, DIM, 0.0);
284
285 // Get the derivatives of the user-defined base state solution
286 // velocity components w.r.t. global coordinates
287 get_base_flow_dudx(time, ipt, interpolated_x, base_flow_dudx);
288
289
290 // MOMENTUM EQUATIONS
291 //------------------
292
293 // Number of master nodes
294 unsigned n_master = 1;
295
296 // Storage for the weight of the shape function
297 double hang_weight = 1.0;
298
299 // Loop over the test functions
300 for (unsigned l = 0; l < n_node; l++)
301 {
302 // Local boolean to indicate whether the node is hanging
303 bool is_node_hanging = node_pt(l)->is_hanging();
304
305 if (is_node_hanging)
306 {
307 hang_info_pt = node_pt(l)->hanging_pt();
308
309 // Read out number of master nodes from hanging data
310 n_master = hang_info_pt->nmaster();
311 }
312 // Otherwise the node is its own master
313 else
314 {
315 n_master = 1;
316 }
317
318 // Loop over the master nodes
319 for (unsigned m = 0; m < n_master; m++)
320 {
321 // Loop over the velocity components
322 for (unsigned i = 0; i < DIM; i++)
323 {
324 // Assemble the residuals
325 // Time dependent term
326 std::complex<double> residual_contribution =
327 -scaled_re_st * eigenvalue * interpolated_u[i] * testf[l] * W;
328 // Pressure term
329 residual_contribution += interpolated_p * dtestfdx(l, i) * W;
330 // Viscous terms
331 for (unsigned k = 0; k < DIM; ++k)
332 {
333 residual_contribution -=
334 visc_ratio *
335 (interpolated_dudx[i][k] + Gamma[i] * interpolated_dudx[k][i]) *
336 dtestfdx(l, k) * W;
337 }
338
339 // Advective terms
340 for (unsigned k = 0; k < DIM; ++k)
341 {
342 residual_contribution -=
343 scaled_re *
344 (base_flow_u[k] * interpolated_dudx[i][k] +
345 interpolated_u[k] * base_flow_dudx(i, k)) *
346 testf[l] * W;
347 }
348
349 // Now separate real and imaginary parts
350
351 if (is_node_hanging)
352 {
353 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
354 u_nodal_index[2 * i]);
355 hang_weight = hang_info_pt->master_weight(m);
356 }
357 // If node is not hanging number or not
358 else
359 {
360 local_eqn = nodal_local_eqn(l, u_nodal_index[2 * i]);
361 hang_weight = 1.0;
362 }
363
364 if (local_eqn >= 0)
365 {
366 residuals[local_eqn] +=
367 residual_contribution.real() * hang_weight;
368 }
369
370
371 if (is_node_hanging)
372 {
373 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
374 u_nodal_index[2 * i + 1]);
375 hang_weight = hang_info_pt->master_weight(m);
376 }
377 // If node is not hanging number or not
378 else
379 {
380 local_eqn = nodal_local_eqn(l, u_nodal_index[2 * i + 1]);
381 hang_weight = 1.0;
382 }
383 if (local_eqn >= 0)
384 {
385 residuals[local_eqn] +=
386 residual_contribution.imag() * hang_weight;
387 }
388
389
390 // CALCULATE THE JACOBIAN
391 /*if(flag)
392 {
393 //Loop over the velocity shape functions again
394 for(unsigned l2=0;l2<n_node;l2++)
395 {
396 //Loop over the velocity components again
397 for(unsigned i2=0;i2<DIM;i2++)
398 {
399 //If at a non-zero degree of freedom add in the entry
400 local_unknown = nodal_local_eqn(l2,u_nodal_index[i2]);
401 if(local_unknown >= 0)
402 {
403 //Add contribution to Elemental Matrix
404 jacobian(local_eqn,local_unknown)
405 -= visc_ratio*Gamma[i]*dpsifdx(l2,i)*dtestfdx(l,i2)*W;
406
407 //Extra component if i2 = i
408 if(i2 == i)
409 {
410 //Loop over velocity components
411 for(unsigned k=0;k<DIM;k++)
412 {
413 jacobian(local_eqn,local_unknown)
414 -= visc_ratio*dpsifdx(l2,k)*dtestfdx(l,k)*W;
415 }
416 }
417
418 //Now add in the inertial terms
419 jacobian(local_eqn,local_unknown)
420 -= scaled_re*psif[l2]*interpolated_dudx(i,i2)*testf[l]*W;
421
422 //Extra component if i2=i
423 if(i2 == i)
424 {
425 //Add the mass matrix term (only diagonal entries)
426 //Note that this is positive because the mass matrix
427 //is taken to the other side of the equation when
428 //formulating the generalised eigenproblem.
429 if(flag==2)
430 {
431 mass_matrix(local_eqn,local_unknown) +=
432 scaled_re_st*psif[l2]*testf[l]*W;
433 }
434
435 //du/dt term
436 jacobian(local_eqn,local_unknown)
437 -= scaled_re_st*
438 node_pt(l2)->time_stepper_pt()->weight(1,0)*
439 psif[l2]*testf[l]*W;
440
441 //Loop over the velocity components
442 for(unsigned k=0;k<DIM;k++)
443 {
444 double tmp=scaled_re*interpolated_u[k];
445 if (!ALE_is_disabled) tmp-=scaled_re_st*mesh_velocity[k];
446 jacobian(local_eqn,local_unknown) -=
447 tmp*dpsifdx(l2,k)*testf[l]*W;
448 }
449 }
450
451 }
452 }
453 }
454
455 //Now loop over pressure shape functions
456 //This is the contribution from pressure gradient
457 for(unsigned l2=0;l2<n_pres;l2++)
458 {
459 //If we are at a non-zero degree of freedom in the entry
460 local_unknown = p_local_eqn(l2);
461 if(local_unknown >= 0)
462 {
463 jacobian(local_eqn,local_unknown)
464 += psip[l2]*dtestfdx(l,i)*W;
465 }
466 }
467 } //End of Jacobian calculation
468
469 }*/ //End of if not boundary condition statement
470
471 } // End of loop over dimension
472 } // End of loop over master nodes
473 } // End of loop over shape functions
474
475
476 // CONTINUITY EQUATION
477 //-------------------
478
479 // Loop over the shape functions
480 for (unsigned l = 0; l < n_pres; l++)
481 {
482 if (pressure_dof_is_hanging[l])
483 {
484 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index[0]);
485 n_master = hang_info_pt->nmaster();
486 }
487 else
488 {
489 n_master = 1;
490 }
491
492 // Loop over the master nodes
493 for (unsigned m = 0; m < n_master; ++m)
494 {
495 // Assemble the residuals
496 std::complex<double> residual_contribution = interpolated_dudx[0][0];
497 for (unsigned k = 1; k < DIM; ++k)
498 {
499 residual_contribution += interpolated_dudx[k][k];
500 }
501
502 if (pressure_dof_is_hanging[l])
503 {
504 local_eqn =
505 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index[0]);
506 hang_weight = hang_info_pt->master_weight(m);
507 }
508 else
509 {
510 local_eqn = this->p_local_eqn(l, 0);
511 hang_weight = 1.0;
512 }
513
514 // If not a boundary conditions
515 if (local_eqn >= 0)
516 {
517 residuals[local_eqn] +=
518 residual_contribution.real() * testp[l] * W * hang_weight;
519 }
520
521 if (pressure_dof_is_hanging[l])
522 {
523 local_eqn =
524 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index[1]);
525 hang_weight = hang_info_pt->master_weight(m);
526 }
527 else
528 {
529 local_eqn = this->p_local_eqn(l, 1);
530 hang_weight = 1.0;
531 }
532
533 // If not a boundary conditions
534 if (local_eqn >= 0)
535 {
536 residuals[local_eqn] +=
537 residual_contribution.imag() * testp[l] * W * hang_weight;
538 }
539
540 } // End of loop over master nodes
541 } // End of loop over l
542
543 // Normalisation condition. Leave this alone because there is
544 // no test function involved.
545 std::complex<double> residual_contribution =
546 interpolated_p_normalisation * interpolated_p;
547 for (unsigned k = 0; k < DIM; ++k)
548 {
549 residual_contribution +=
550 interpolated_u_normalisation[k] * interpolated_u[k];
551 }
552
553 local_eqn = this->eigenvalue_local_eqn(0);
554 if (local_eqn >= 0)
555 {
556 residuals[local_eqn] += residual_contribution.real() * W;
557 }
558
559 local_eqn = this->eigenvalue_local_eqn(1);
560 if (local_eqn >= 0)
561 {
562 residuals[local_eqn] += residual_contribution.imag() * W;
563 }
564
565 } // End of loop over the integration points
566
567 } // End of fill_in_generic_residual_contribution_linearised_nst
568
569} // End of namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
const double & viscosity_ratio() const
Viscosity ratio for element: element's viscosity relative to the viscosity used in the definition of ...
virtual unsigned u_index_linearised_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
virtual unsigned npres_linearised_nst() const =0
Return the number of pressure degrees of freedom associated with a single pressure component in the e...
const double & density_ratio() const
Density ratio for element: element's density relative to the viscosity used in the definition of the ...
virtual int p_local_eqn(const unsigned &n, const unsigned &i)=0
Access function for the local equation number information for the i-th component of the pressure....
virtual void pshape_linearised_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
virtual void get_base_flow_u(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
Calculate the velocity components of the base flow solution at a given time and Eulerian position.
virtual double p_linearised_nst(const unsigned &n_p, const unsigned &i) const =0
Return the i-th pressure value at local pressure "node" n_p. Uses suitably interpolated value for han...
virtual double dshape_and_dtest_eulerian_at_knot_linearised_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and their derivatives w.r.t. global coordinates at the ipt-th integration...
virtual void get_base_flow_dudx(const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
Calculate the derivatives of the velocity components of the base flow solution w.r....
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
virtual int p_index_linearised_nst(const unsigned &i) const
Which nodal value represents the pressure?
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
void fill_in_generic_residual_contribution_linearised_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to the elemental residual vector and/or Jacobian matrix flag=1: compute bo...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
double & time()
Return current value of continous time.
Definition: timesteppers.h:332
//////////////////////////////////////////////////////////////////// ////////////////////////////////...