axisymmetric_advection_navier_stokes_elements.h
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// Version 1.0; svn revision $LastChangedRevision: 1162 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-04-18 13:27:54 +0100 (Mon, 18 Apr 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef AXISYMMETRIC_ADVECTION_NAVIER_STOKES_ELEMENTS_HEADER
31 #define AXISYMMETRIC_ADVECTION_NAVIER_STOKES_ELEMENTS_HEADER
32 
33 //Oomph-lib headers, we require the generic, advection-diffusion-reaction,
34 //navier-stokes elements and fluid interface elements
35 #include "generic.h"
36 #include "axisym_advection_diffusion.h"
37 #include "axisym_navier_stokes.h"
38 
39 namespace oomph
40 {
41 
42 //======================class definition==============================
43 ///A class that solves the Boussinesq approximation of the Navier--Stokes
44 ///and energy equations by coupling two pre-existing classes.
45 ///The QAdvectionDiffusionReactionElement with
46 ///bi-quadratic interpolation for the
47 ///scalar variables (temperature and concentration) and
48 ///QCrouzeixRaviartElement which solves the Navier--Stokes equations
49 ///using bi-quadratic interpolation for the velocities and a discontinuous
50 ///bi-linear interpolation for the pressure. Note that we are free to
51 ///choose the order in which we store the variables at the nodes. In this
52 ///case we choose to store the variables in the order fluid velocities
53 ///followed by temperature. We must, therefore, overload the function
54 ///AdvectionDiffusionReactionEquations<2,DIM>::u_index_adv_diff_react()
55 ///to indicate that
56 ///the temperature is stored at the DIM-th position not the 0-th. We do not
57 ///need to overload the corresponding function in the
58 ///NavierStokesEquations<DIM> class because the velocities are stored
59 ///first.
60 //=========================================================================
62  public virtual QAxisymAdvectionDiffusionElement<3>,
63  public virtual AxisymmetricQCrouzeixRaviartElement
64 {
65 
66 public:
67 
68  /// \short Constructor: call the underlying constructors and
69  /// initialise the pointer to the Rayleigh number to point
70  /// to the default value of 0.0.
72  QAxisymAdvectionDiffusionElement<3>(),
73  AxisymmetricQCrouzeixRaviartElement() {}
74 
75 
76  ///\short The required number of values stored at the nodes is the sum of the
77  ///required values of the two single-physics elements. Note that this step is
78  ///generic for any multi-physics element of this type.
79  unsigned required_nvalue(const unsigned &n) const
80  {return (QAxisymAdvectionDiffusionElement<3>::required_nvalue(n) +
81  AxisymmetricQCrouzeixRaviartElement::required_nvalue(n));}
82 
83  /// Final override for disable ALE
84  void disable_ALE()
85  {
86  //Disable ALE in both sets of equations
87  AxisymmetricNavierStokesEquations::disable_ALE();
88  AxisymAdvectionDiffusionEquations::enable_ALE();
89  }
90 
91  /// Final override for enable ALE
92  void enable_ALE()
93  {
94  //Enable ALE in both sets of equations
95  AxisymmetricNavierStokesEquations::enable_ALE();
96  AxisymAdvectionDiffusionEquations::enable_ALE();
97  }
98 
99 
100  /// \short Number of scalars/fields output by this element. Reimplements
101  /// broken virtual function in base class.
102  unsigned nscalar_paraview() const
103  {
104  return AxisymmetricNavierStokesEquations::nscalar_paraview()
105  + AxisymAdvectionDiffusionEquations::nscalar_paraview();
106  }
107 
108  /// \short Write values of the i-th scalar field at the plot points. Needs
109  /// to be implemented for each new specific element type.
110  void scalar_value_paraview(std::ofstream& file_out,
111  const unsigned& i,
112  const unsigned& nplot) const
113  {
114  AxisymmetricNavierStokesEquations::scalar_value_paraview(file_out,i,nplot);
115  AxisymAdvectionDiffusionEquations::scalar_value_paraview(file_out,i,nplot);
116  }
117 
118 
119  std::string scalar_name_paraview(const unsigned& i) const
120  {
121  return AxisymmetricNavierStokesEquations::scalar_name_paraview(i);
122  //AxisymAdvectionDiffusionEquations::scalar_name_paraview(i);
123  }
124 
125 
126 
127  /// Overload the standard output function with the broken default
128  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
129 
130  /// \short Output function:
131  /// Output x, y, u, v, p, theta at Nplot^DIM plot points
132  // Start of output function
133  void output(std::ostream &outfile, const unsigned &nplot)
134  {
135  //vector of local coordinates
136  Vector<double> s(2);
137 
138  // Tecplot header info
139  outfile << this->tecplot_zone_string(nplot);
140 
141  // Loop over plot points
142  unsigned num_plot_points=this->nplot_points(nplot);
143  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
144  {
145  // Get local coordinates of plot point
146  this->get_s_plot(iplot,nplot,s);
147 
148  // Output the position of the plot point
149  for(unsigned i=0;i<2;i++)
150  {outfile << this->interpolated_x(s,i) << " ";}
151 
152  // Output the fluid velocities at the plot point
153  for(unsigned i=0;i<3;i++)
154  {outfile << this->interpolated_u_axi_nst(s,i) << " ";}
155 
156  // Output the fluid pressure at the plot point
157  outfile << this->interpolated_p_axi_nst(s) << " ";
158 
159  //Output the solute concentration
160  outfile << this->interpolated_u_axi_adv_diff(s) << " ";
161 
162  outfile << "\n";
163  }
164  outfile << std::endl;
165 
166  // Write tecplot footer (e.g. FE connectivity lists)
167  this->write_tecplot_zone_footer(outfile,nplot);
168  } //End of output function
169 
170 
171  /// \short C-style output function: Broken default
172  void output(FILE* file_pt)
173  {FiniteElement::output(file_pt);}
174 
175  /// \short C-style output function: Broken default
176  void output(FILE* file_pt, const unsigned &n_plot)
177  {FiniteElement::output(file_pt,n_plot);}
178 
179  /// \short Output function for an exact solution: Broken default
180  void output_fct(std::ostream &outfile, const unsigned &Nplot,
181  FiniteElement::SteadyExactSolutionFctPt
182  exact_soln_pt)
183  {FiniteElement::output_fct(outfile,Nplot,exact_soln_pt);}
184 
185 
186  /// \short Output function for a time-dependent exact solution:
187  /// Broken default.
188  void output_fct(std::ostream &outfile, const unsigned &Nplot,
189  const double& time,
190  FiniteElement::UnsteadyExactSolutionFctPt
191  exact_soln_pt)
192  {
193  FiniteElement::
194  output_fct(outfile,Nplot,time,exact_soln_pt);
195  }
196 
197  ///\short Overload the index at which the temperature and solute
198  ///concentration variables are stored.
199  // We choose to store them after the fluid velocities.
200  inline unsigned u_index_axi_adv_diff() const {return 3;}
201 
202  /// \short Validate against exact solution at given time
203  /// Solution is provided via function pointer.
204  /// Plot at a given number of plot points and compute L2 error
205  /// and L2 norm of velocity solution over element
206  /// Call the broken default
207  void compute_error(std::ostream &outfile,
208  FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
209  const double& time,
210  double& error, double& norm)
211  {FiniteElement::compute_error(outfile,exact_soln_pt,
212  time,error,norm);}
213 
214  /// \short Validate against exact solution.
215  /// Solution is provided via function pointer.
216  /// Plot at a given number of plot points and compute L2 error
217  /// and L2 norm of velocity solution over element
218  /// Call the broken default
219  void compute_error(std::ostream &outfile,
220  FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
221  double& error, double& norm)
222  {FiniteElement::compute_error(outfile,exact_soln_pt,error,norm);}
223 
224  /// \short Overload the wind function in the advection-diffusion equations.
225  /// This provides the coupling from the Navier--Stokes equations to the
226  /// advection-diffusion equations because the wind is the fluid velocity.
227  void get_wind_axi_adv_diff(const unsigned& ipt,
228  const Vector<double> &s, const Vector<double>& x,
229  Vector<double>& wind) const
230  {
231  //The wind function is simply the velocity at the points
232  this->interpolated_u_axi_nst(s,wind);
233  }
234 
235  /// \short Calculate the element's contribution to the residual vector.
236  /// Recall that fill_in_* functions MUST NOT initialise the entries
237  /// in the vector to zero. This allows us to call the
238  /// fill_in_* functions of the constituent single-physics elements
239  /// sequentially, without wiping out any previously computed entries.
240  void fill_in_contribution_to_residuals(Vector<double> &residuals)
241  {
242  //Fill in the residuals of the Navier-Stokes equations
243  AxisymmetricNavierStokesEquations::fill_in_contribution_to_residuals(residuals);
244 
245  //Fill in the residuals of the advection-diffusion eqautions
246  AxisymAdvectionDiffusionEquations::fill_in_contribution_to_residuals(residuals);
247  }
248 
249 
250 #ifdef USE_FD_JACOBIAN_FOR_ADVECTION_DIFFUSION_NAVIER_STOKES_ELEMENT
251 
252 
253  ///\short Compute the element's residual vector and the Jacobian matrix.
254  /// Jacobian is computed by finite-differencing.
255  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
256  DenseMatrix<double> &jacobian)
257  {
258  // This function computes the Jacobian by finite-differencing
259  FiniteElement::fill_in_contribution_to_jacobian(residuals,jacobian);
260  }
261 
262 #else
263 
264  ///\short Helper function to get the off-diagonal blocks of the Jacobian
265  ///matrix by finite differences
266  void fill_in_off_diagonal_jacobian_blocks_by_fd(Vector<double> &residuals,
267  DenseMatrix<double> &jacobian)
268  {
269  //Local storage for the index in the nodes at which the
270  //Navier-Stokes velocities are stored (we know that this should be 0,1,2)
271  unsigned u_nodal_nst[3];
272  for(unsigned i=0;i<3;i++)
273  {u_nodal_nst[i] = this->u_index_nst(i);}
274 
275  //Local storage for the index at which the temperature and
276  //solute are stored
277  unsigned C_nodal_adv_diff = this->u_index_axi_adv_diff();
278 
279 
280  //Find the total number of unknowns in the elements
281  unsigned n_dof = this->ndof();
282 
283  //Temporary storage for residuals
284  Vector<double> newres(n_dof);
285 
286  //Storage for local unknown and local equation
287  int local_unknown =0, local_eqn = 0;
288 
289  //Set the finite difference step
290  double fd_step = FiniteElement::Default_fd_jacobian_step;
291 
292  //Find the number of nodes
293  unsigned n_node = this->nnode();
294 
295  //Calculate the contribution of the Navier--Stokes velocities to the
296  //advection-diffusion equations
297 
298  //Loop over the nodes
299  for(unsigned n=0;n<n_node;n++)
300  {
301  //There are 3 values of the velocities
302  for(unsigned i=0;i<3;i++)
303  {
304  //Get the local velocity equation number
305  local_unknown = this->nodal_local_eqn(n,u_nodal_nst[i]);
306 
307  //If it's not pinned
308  if(local_unknown >= 0)
309  {
310  //Get a pointer to the velocity value
311  double *value_pt = this->node_pt(n)->value_pt(u_nodal_nst[i]);
312 
313  //Save the old value
314  double old_var = *value_pt;
315 
316  //Increment the value
317  *value_pt += fd_step;
318 
319  //Get the altered advection-diffusion residuals.
320  //Do this using fill_in because get_residuals has never been
321  //overloaded, and will actually compute all residuals which
322  //is slightly inefficient.
323  for(unsigned m=0;m<n_dof;m++) {newres[m] = 0.0;}
324  AxisymAdvectionDiffusionEquations::
325  fill_in_contribution_to_residuals(newres);
326 
327  //Now fill in the Advection-Diffusion-Reaction sub-block
328  //of the jacobian
329  for(unsigned m=0;m<n_node;m++)
330  {
331  //Local equation for temperature or solute
332  local_eqn = this->nodal_local_eqn(m,C_nodal_adv_diff);
333 
334  //If it's not a boundary condition
335  if(local_eqn >= 0)
336  {
337  double sum = (newres[local_eqn] - residuals[local_eqn])/fd_step;
338  jacobian(local_eqn,local_unknown) = sum;
339  }
340  }
341 
342  //Reset the nodal data
343  *value_pt = old_var;
344  }
345  }
346 
347  } //End of loop over nodes
348 
349 }
350 
351  ///\short Compute the element's residual Vector and the Jacobian matrix.
352  /// Use finite-differencing only for the off-diagonal blocks.
353  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
354  DenseMatrix<double> &jacobian)
355  {
356 
357  //Calculate the Navier-Stokes contributions (diagonal block and residuals)
358  AxisymmetricNavierStokesEquations::
359  fill_in_contribution_to_jacobian(residuals,jacobian);
360 
361  //Calculate the advection-diffusion contributions
362  //(diagonal block and residuals)
363  AxisymAdvectionDiffusionEquations::
364  fill_in_contribution_to_jacobian(residuals,jacobian);
365 
366  //Add in the off-diagonal blocks
367  fill_in_off_diagonal_jacobian_blocks_by_fd(residuals,jacobian);
368  } //End of jacobian calculation
369 
370 #endif
371 
372 
373  /// Add the element's contribution to its residuals vector,
374  /// jacobian matrix and mass matrix
376  Vector<double> &residuals, DenseMatrix<double> &jacobian,
377  DenseMatrix<double> &mass_matrix)
378  {
379  AxisymmetricNavierStokesEquations::
380  fill_in_contribution_to_jacobian_and_mass_matrix(
381  residuals,jacobian,mass_matrix);
382 
383  AxisymAdvectionDiffusionEquations::
384  fill_in_contribution_to_jacobian_and_mass_matrix(
385  residuals,jacobian,mass_matrix);
386 
387  fill_in_off_diagonal_jacobian_blocks_by_fd(residuals,jacobian);
388  }
389 
390 
391 };
392 
393 
394 //=======================================================================
395 /// Face geometry of the 2D DoubleBuoyant Crouzeix_Raviart elements
396 //=======================================================================
397 template<>
399 public virtual QElement<1,3>
400 {
401  public:
402  FaceGeometry() : QElement<1,3>() {}
403 };
404 
405 //=======================================================================
406 /// Face geometry of the 2D DoubleBuoyant Crouzeix_Raviart elements
407 //=======================================================================
408 template<>
409 class FaceGeometry<FaceGeometry<AxisymmetricQAdvectionCrouzeixRaviartElement> >:
410 public virtual PointElement
411 {
412  public:
413  FaceGeometry() : PointElement() {}
414 };
415 
416 } //End of the oomph namespace
417 
418 #endif
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: Broken default.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: Output x, y, u, v, p, theta at Nplot^DIM plot points.
void fill_in_off_diagonal_jacobian_blocks_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Helper function to get the off-diagonal blocks of the Jacobian matrix by finite differences.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Calculate the element's contribution to the residual vector. Recall that fill_in_* functions MUST NOT...
void output(FILE *file_pt)
C-style output function: Broken default.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
AxisymmetricQAdvectionCrouzeixRaviartElement()
Constructor: call the underlying constructors and initialise the pointer to the Rayleigh number to po...
void output_fct(std::ostream &outfile, const unsigned &Nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: Broken default.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Validate against exact solution. Solution is provided via function pointer. Plot at a given number of...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual vector and the Jacobian matrix. Jacobian is computed by finite-differe...
unsigned u_index_axi_adv_diff() const
Overload the index at which the temperature and solute concentration variables are stored...
void output(std::ostream &outfile)
Overload the standard output function with the broken default.
void output_fct(std::ostream &outfile, const unsigned &Nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution: Broken default.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element Call the broken default.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
unsigned required_nvalue(const unsigned &n) const
The required number of values stored at the nodes is the sum of the required values of the two single...
void get_wind_axi_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Overload the wind function in the advection-diffusion equations. This provides the coupling from the ...