spherical_advection_diffusion_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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 //Non-inline functions for Advection Diffusion elements in a spherical
31 //coordinate system
33 
34 namespace oomph
35 {
36 
37 /// 2D Steady Axisymmetric Advection Diffusion elements
38 
39 /// Default value for Peclet number
42 
43 
44 //======================================================================
45 /// \short Compute element residual Vector and/or element Jacobian matrix
46 ///
47 /// flag=1: compute both
48 /// flag=0: compute only residual Vector
49 ///
50 /// Pure version without hanging nodes
51 //======================================================================
54  Vector<double> &residuals, DenseMatrix<double> &jacobian,
55  DenseMatrix<double> &mass_matrix,
56  unsigned flag)
57 {
58  //Find out how many nodes there are
59  const unsigned n_node = nnode();
60 
61  //Get the nodal index at which the unknown is stored
62  const unsigned u_nodal_index = u_index_spherical_adv_diff();
63 
64  //Set up memory for the shape and test functions
65  Shape psi(n_node), test(n_node);
66  DShape dpsidx(n_node,2), dtestdx(n_node,2);
67 
68  //Set the value of n_intpt
69  const unsigned n_intpt = integral_pt()->nweight();
70 
71  //Set the Vector to hold local coordinates
72  Vector<double> s(2);
73 
74  //Get Physical Variables from Element
75  const double scaled_peclet = pe();
76 
77  //Get peclet number multiplied by Strouhal number
78  const double scaled_peclet_st = pe_st();
79 
80  //Integers used to store the local equation number and local unknown
81  //indices for the residuals and jacobians
82  int local_eqn=0, local_unknown=0;
83 
84  //Loop over the integration points
85  for(unsigned ipt=0;ipt<n_intpt;ipt++)
86  {
87 
88  //Assign values of s
89  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
90 
91  //Get the integral weight
92  double w = integral_pt()->weight(ipt);
93 
94  //Call the derivatives of the shape and test functions
95  double J =
97  ipt,psi,dpsidx,test,dtestdx);
98 
99  //Premultiply the weights and the Jacobian
100  double W = w*J;
101 
102  //Calculate local values of the solution and its derivatives
103  //Allocate
104  double interpolated_u = 0.0;
105  double dudt = 0.0;
107  Vector<double> interpolated_dudx(2,0.0);
108 
109  //Calculate function value and derivatives:
110  //-----------------------------------------
111  //Loop over nodes
112  for(unsigned l=0;l<n_node;l++)
113  {
114  //Get the value at the node
115  double u_value = raw_nodal_value(l,u_nodal_index);
116  interpolated_u += u_value*psi(l);
117  dudt += du_dt_spherical_adv_diff(l)*psi(l);
118  //Loop over the two coordinates directions
119  for(unsigned j=0;j<2;j++)
120  {
121  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
122  interpolated_dudx[j] += u_value*dpsidx(l,j);
123  }
124  }
125 
126  //Get source function
127  //-------------------
128  double source;
129  get_source_spherical_adv_diff(ipt,interpolated_x,source);
130 
131 
132  //Get wind
133  //--------
134  Vector<double> wind(3);
135  get_wind_spherical_adv_diff(ipt,s,interpolated_x,wind);
136 
137  //r is the first position component
138  double r = interpolated_x[0];
139  //theta is the second position component
140  double sin_th = sin(interpolated_x[1]);
141  //dS is the area weighting
142  double dS = r*r*sin_th;
143 
144  //TEMPERATURE EQUATION (Neglected viscous dissipation)
145  //Assemble residuals and Jacobian
146  //--------------------------------
147 
148  // Loop over the test functions
149  for(unsigned l=0;l<n_node;l++)
150  {
151  //Set the local equation number
152  local_eqn = nodal_local_eqn(l,u_nodal_index);
153 
154  /*IF it's not a boundary condition*/
155  if(local_eqn >= 0)
156  {
157  //Add body force/source term
158  residuals[local_eqn] -= (scaled_peclet_st*dudt + source)*dS*test(l)*W;
159 
160  //The Advection Diffusion bit itself
161  residuals[local_eqn] -=
162  //radial terms
163  (dS*interpolated_dudx[0]*
164  (scaled_peclet*wind[0]*test(l) + dtestdx(l,0)) +
165  //azimuthal terms
166  (sin_th*interpolated_dudx[1]*
167  (r*scaled_peclet*wind[1]*test(l) + dtestdx(l,1))))*W;
168 
169 
170  // Calculate the jacobian
171  //-----------------------
172  if(flag)
173  {
174  //Loop over the velocity shape functions again
175  for(unsigned l2=0;l2<n_node;l2++)
176  {
177  //Set the number of the unknown
178  local_unknown = nodal_local_eqn(l2,u_nodal_index);
179 
180  //If at a non-zero degree of freedom add in the entry
181  if(local_unknown >= 0)
182  {
183  //Mass matrix term
184  jacobian(local_eqn,local_unknown)
185  -= scaled_peclet_st*test(l)*psi(l2)*
186  node_pt(l2)->time_stepper_pt()->weight(1,0)*dS*W;
187 
188  //add the mass matrix term
189  if(flag==2)
190  {
191  mass_matrix(local_eqn,local_unknown)
192  += scaled_peclet_st*test(l)*psi(l2)*dS*W;
193  }
194 
195  //Assemble Jacobian term
196  jacobian(local_eqn,local_unknown) -=
197  //radial terms
198  (dS*dpsidx(l2,0)*
199  (scaled_peclet*wind[0]*test(l) + dtestdx(l,0)) +
200  //azimuthal terms
201  (sin_th*dpsidx(l2,1)*
202  (r*scaled_peclet*wind[1]*test(l) + dtestdx(l,1))))*W;
203  }
204  }
205  }
206  }
207  }
208 
209 
210  } // End of loop over integration points
211 }
212 
213 
214 
215 
216 //======================================================================
217 /// Self-test: Return 0 for OK
218 //======================================================================
219 //template <unsigned DIM>
221 {
222 
223  bool passed = true;
224 
225  // Check lower-level stuff
226  if (FiniteElement::self_test()!=0)
227  {
228  passed = false;
229  }
230 
231  // Return verdict
232  if (passed)
233  {
234  return 0;
235  }
236  else
237  {
238  return 1;
239  }
240 
241 }
242 
243 
244 
245 //======================================================================
246 /// \short Output function:
247 ///
248 /// r,z,u,w_r,w_z
249 ///
250 /// nplot points in each coordinate direction
251 //======================================================================
253  const unsigned &nplot)
254 {
255  //Vector of local coordinates
256  Vector<double> s(2);
257 
258  //Tecplot header info
259  outfile << tecplot_zone_string(nplot);
260 
261  //Loop over plot points
262  unsigned num_plot_points = nplot_points(nplot);
263  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
264  {
265  //Get local coordinates of plot point
266  get_s_plot(iplot,nplot,s);
267 
268  //Get Eulerian coordinate of plot point
269  Vector<double> x(2);
270  interpolated_x(s,x);
271 
272  for(unsigned i=0;i<2;i++)
273  {
274  outfile << x[i] << " ";
275  }
276 
277  //Convert to Cartesians
278  outfile << x[0]*sin(x[1]) << " " << x[0]*cos(x[1]) << " ";
279  //Output concentration
280  outfile << this->interpolated_u_spherical_adv_diff(s) << " ";
281 
282  //Get the wind
283  Vector<double> wind(2);
284  //Dummy ipt argument needed... ?
285  unsigned ipt = 0;
286  get_wind_spherical_adv_diff(ipt,s,x,wind);
287  for(unsigned i=0;i<2;i++)
288  {
289  outfile << wind[i] << " ";
290  }
291  outfile << std::endl;
292 
293  }
294 
295  //Write tecplot footer (e.g. FE connectivity lists)
296  write_tecplot_zone_footer(outfile,nplot);
297 
298 }
299 
300 
301 //======================================================================
302 /// C-style output function:
303 ///
304 /// r,z,u
305 ///
306 /// nplot points in each coordinate direction
307 //======================================================================
308 //template <unsigned DIM>
310  const unsigned &nplot)
311 {
312  //Vector of local coordinates
313  Vector<double> s(2);
314 
315  //Tecplot header info
316  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
317 
318  //Loop over plot points
319  unsigned num_plot_points = nplot_points(nplot);
320  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
321  {
322 
323  //Get local coordinates of plot point
324  get_s_plot(iplot,nplot,s);
325 
326  for(unsigned i=0;i<2;i++)
327  {
328  fprintf(file_pt,"%g ",interpolated_x(s,i));
329 
330  }
331  //Write Cartesian coordinates
332  double r = interpolated_x(s,0); double theta = interpolated_x(s,1);
333  fprintf(file_pt,"%g %g ",r*sin(theta),r*cos(theta));
334 
335  fprintf(file_pt,"%g \n",interpolated_u_spherical_adv_diff(s));
336  }
337 
338  //Write tecplot footer (e.g. FE connectivity lists)
339  write_tecplot_zone_footer(file_pt,nplot);
340 
341 }
342 
343 //======================================================================
344  /// \short Output exact solution
345  ///
346  /// Solution is provided via function pointer.
347  /// Plot at a given number of plot points.
348  ///
349  /// r,z,u_exact
350 //======================================================================
351 //template <unsigned DIM>
353  const unsigned &nplot,
355  {
356 
357  //Vector of local coordinates
358  Vector<double> s(2);
359 
360  //Vector for coordintes
361  Vector<double> x(2);
362 
363  //Tecplot header info
364  outfile << tecplot_zone_string(nplot);
365 
366  //Exact solution Vector (here a scalar)
367  Vector<double> exact_soln(1);
368 
369  //Loop over plot points
370  unsigned num_plot_points=nplot_points(nplot);
371  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
372  {
373 
374  //Get local coordinates of plot point
375  get_s_plot(iplot,nplot,s);
376 
377  //Get x position as Vector
378  interpolated_x(s,x);
379 
380  //Get exact solution at this point
381  (*exact_soln_pt)(x,exact_soln);
382 
383  //Output x,y,...,u_exact
384  for(unsigned i=0;i<2;i++)
385  {
386  outfile << x[i] << " ";
387  }
388  outfile << exact_soln[0] << std::endl;
389  }
390 
391  //Write tecplot footer (e.g. FE connectivity lists)
392  write_tecplot_zone_footer(outfile,nplot);
393 
394  }
395 
396 
397 //======================================================================
398  /// \short Validate against exact solution
399  ///
400  /// Solution is provided via function pointer.
401  /// Plot error at a given number of plot points.
402  ///
403 //======================================================================
404 //template <unsigned DIM>
406  std::ostream &outfile,
408  double& error,
409  double& norm)
410 {
411 
412  //Initialise
413  error=0.0;
414  norm=0.0;
415 
416  //Vector of local coordinates
417  Vector<double> s(2);
418 
419  //Vector for coordintes
420  Vector<double> x(2);
421 
422  //Find out how many nodes there are in the element
423  unsigned n_node = nnode();
424 
425  Shape psi(n_node);
426 
427  //Set the value of n_intpt
428  unsigned n_intpt = integral_pt()->nweight();
429 
430  //Tecplot header info
431  outfile << "ZONE" << std::endl;
432 
433  //Exact solution Vector (here a scalar)
434  Vector<double> exact_soln(1);
435 
436  //Loop over the integration points
437  for(unsigned ipt=0;ipt<n_intpt;ipt++)
438  {
439 
440  //Assign values of s
441  for(unsigned i=0;i<2;i++)
442  {
443  s[i] = integral_pt()->knot(ipt,i);
444  }
445 
446  //Get the integral weight
447  double w = integral_pt()->weight(ipt);
448 
449  //Get jacobian of mapping
450  double J=J_eulerian(s);
451 
452  //Premultiply the weights and the Jacobian
453  double W = w*J;
454 
455  //Get x position as Vector
456  interpolated_x(s,x);
457 
458  //Get FE function value
459  double u_fe=interpolated_u_spherical_adv_diff(s);
460 
461  //Get exact solution at this point
462  (*exact_soln_pt)(x,exact_soln);
463 
464  //Output x,y,...,error
465  for(unsigned i=0;i<2;i++)
466  {
467  outfile << x[i] << " ";
468  }
469  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
470 
471  //Add to error and norm
472  norm+=exact_soln[0]*exact_soln[0]*W;
473  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
474 
475  }
476 
477 }
478 
479 
480 //======================================================================
481 // Set the data for the number of Variables at each node
482 //======================================================================
483 template<unsigned NNODE_1D>
485 Initial_Nvalue = 1;
486 
487 //====================================================================
488 // Force build of templates
489 //====================================================================
490 
494 
495 
496 }
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
void output(std::ostream &outfile)
Output with default number of plot points.
const double & pe_st() const
Peclet number multiplied by Strouhal number.
cstr elem_len * i
Definition: cfortran.h:607
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
double interpolated_u_spherical_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
virtual void fill_in_generic_residual_contribution_spherical_adv_diff(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add the element's contribution to its residual vector only (if flag=and/or element Jacobian matrix...
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation. NOTE: Moved to cc file because of a possible compiler bug in gcc (yes, really!). The move to the cc file avoids inlining which appears to cause problems (only) when compiled with gcc and -O3; offensive "illegal read" is in optimised-out section of code and data that is allegedly illegal is readily readable (by other means) just before this function is called so I can't really see how we could possibly be responsible for this...
Definition: elements.cc:1657
static double Default_peclet_number
Static default value for the Peclet number.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
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:3841
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
double du_dt_spherical_adv_diff(const unsigned &n) const
du/dt at local node n.
static char t char * s
Definition: cfortran.h:572
virtual void get_source_spherical_adv_diff(const unsigned &ipt, const Vector< double > &x, double &source) const
Get source term at (Eulerian) position x. This function is virtual to allow overloading in multi-phys...
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,z,u_exact at nplot^2 plot points.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
virtual double dshape_and_dtest_eulerian_at_knot_spherical_adv_diff(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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:1383
QSphericalAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Axisymmetric Advec...
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2446
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
virtual void get_wind_spherical_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Get wind at (Eulerian) position x and/or local coordinate s. This function is virtual to allow overlo...
virtual unsigned u_index_spherical_adv_diff() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...