advection_diffusion_reaction_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
32 
33 namespace oomph
34 {
35 
36 ///2D Advection Diffusion elements
37 
38 
39 /// Default value for Peclet number
40 template<unsigned NREAGENT, unsigned DIM>
42 Default_dimensionless_number(NREAGENT,1.0);
43 
44 //======================================================================
45 /// \short Compute element residual Vector and/or element Jacobian matrix
46 /// and/or the mass matrix
47 ///
48 /// flag=2: compute Jacobian, residuals and mass matrix
49 /// flag=1: compute Jacobian and residuals
50 /// flag=0: compute only residual Vector
51 ///
52 /// Pure version without hanging nodes
53 //======================================================================
54 template <unsigned NREAGENT, unsigned DIM>
57  Vector<double> &residuals,
58  DenseMatrix<double> &jacobian,
60  &mass_matrix,
61  unsigned flag)
62 {
63  //Find out how many nodes there are
64  const unsigned n_node = nnode();
65 
66  //Get the nodal index at which the unknown is stored
67  unsigned c_nodal_index[NREAGENT];
68  for(unsigned r=0;r<NREAGENT;r++)
69  {c_nodal_index[r]= c_index_adv_diff_react(r);}
70 
71  //Set up memory for the shape and test functions
72  Shape psi(n_node), test(n_node);
73  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
74 
75  //Set the value of n_intpt
76  unsigned n_intpt = integral_pt()->nweight();
77 
78  //Set the Vector to hold local coordinates
79  Vector<double> s(DIM);
80 
81  //Get diffusion coefficients
82  Vector<double> D = diff();
83 
84  //Get the timescales
85  Vector<double> T = tau();
86 
87  //Integers used to store the local equation number and local unknown
88  //indices for the residuals and jacobians
89  int local_eqn=0, local_unknown=0;
90 
91  //Loop over the integration points
92  for(unsigned ipt=0;ipt<n_intpt;ipt++)
93  {
94 
95  //Assign values of s
96  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
97 
98  //Get the integral weight
99  double w = integral_pt()->weight(ipt);
100 
101  //Call the derivatives of the shape and test functions
102  double J =
103  dshape_and_dtest_eulerian_at_knot_adv_diff_react(
104  ipt,psi,dpsidx,test,dtestdx);
105 
106  //Premultiply the weights and the Jacobian
107  double W = w*J;
108 
109  //Calculate local values of the solution and its derivatives
110  //Allocate
111  Vector<double> interpolated_c(NREAGENT,0.0);
112  Vector<double> dcdt(NREAGENT,0.0);
113  Vector<double> interpolated_x(DIM,0.0);
114  DenseMatrix<double> interpolated_dcdx(NREAGENT,DIM,0.0);
115  Vector<double> mesh_velocity(DIM,0.0);
116 
117 
118  //Calculate function value and derivatives:
119  // Loop over nodes
120  for(unsigned l=0;l<n_node;l++)
121  {
122  // Loop over directions to calculate the position
123  for(unsigned j=0;j<DIM;j++)
124  {
125  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
126  }
127 
128  //Loop over the unknown reagents
129  for(unsigned r=0;r<NREAGENT;r++)
130  {
131  //Get the value at the node
132  const double c_value = raw_nodal_value(l,c_nodal_index[r]);
133 
134  //Calculate the interpolated value
135  interpolated_c[r] += c_value*psi(l);
136  dcdt[r] += dc_dt_adv_diff_react(l,r)*psi(l);
137 
138  // Loop over directions to calculate the derivatie
139  for(unsigned j=0;j<DIM;j++)
140  {
141  interpolated_dcdx(r,j) += c_value*dpsidx(l,j);
142  }
143  }
144  }
145 
146  // Mesh velocity?
147  if (!ALE_is_disabled)
148  {
149  for(unsigned l=0;l<n_node;l++)
150  {
151  for(unsigned j=0;j<DIM;j++)
152  {
153  mesh_velocity[j] += raw_dnodal_position_dt(l,j)*psi(l);
154  }
155  }
156  }
157 
158 
159  //Get source function
160  Vector<double> source(NREAGENT);
161  get_source_adv_diff_react(ipt,interpolated_x,source);
162 
163 
164  //Get wind
165  Vector<double> wind(DIM);
166  get_wind_adv_diff_react(ipt,s,interpolated_x,wind);
167 
168  //Get reaction terms
169  Vector<double> R(NREAGENT);
170  get_reaction_adv_diff_react(ipt,interpolated_c,R);
171 
172  //If we are getting the jacobian the get the derivative terms
173  DenseMatrix<double> dRdC(NREAGENT);
174  if(flag)
175  {
176  get_reaction_deriv_adv_diff_react(ipt,interpolated_c,dRdC);
177  }
178 
179 
180  // Assemble residuals and Jacobian
181  //--------------------------------
182 
183  // Loop over the test functions
184  for(unsigned l=0;l<n_node;l++)
185  {
186  //Loop over the reagents
187  for(unsigned r=0;r<NREAGENT;r++)
188  {
189  //Set the local equation number
190  local_eqn = nodal_local_eqn(l,c_nodal_index[r]);
191 
192  /*IF it's not a boundary condition*/
193  if(local_eqn >= 0)
194  {
195  // Add body force/source/reaction term and time derivative
196  residuals[local_eqn] -=
197  (T[r]*dcdt[r] + source[r] + R[r])*test(l)*W;
198 
199  // The Advection Diffusion bit itself
200  for(unsigned k=0;k<DIM;k++)
201  {
202  //Terms that multiply the test function
203  double tmp = wind[k];
204  //If the mesh is moving need to subtract the mesh velocity
205  if(!ALE_is_disabled) {tmp -= T[r]*mesh_velocity[k];}
206  //Now construct the contribution to the residuals
207  residuals[local_eqn] -=
208  interpolated_dcdx(r,k)*(tmp*test(l) + D[r]*dtestdx(l,k))*W;
209  }
210 
211  // Calculate the jacobian
212  //-----------------------
213  if(flag)
214  {
215  //Loop over the velocity shape functions again
216  for(unsigned l2=0;l2<n_node;l2++)
217  {
218  //Loop over the reagents again
219  for(unsigned r2=0;r2<NREAGENT;r2++)
220  {
221  //Set the number of the unknown
222  local_unknown = nodal_local_eqn(l2,c_nodal_index[r2]);
223 
224  //If at a non-zero degree of freedom add in the entry
225  if(local_unknown >= 0)
226  {
227  //Diagonal terms (i.e. the basic equations)
228  if(r2==r)
229  {
230  //Mass matrix term
231  jacobian(local_eqn,local_unknown)
232  -= T[r]*test(l)*psi(l2)*
233  node_pt(l2)->time_stepper_pt()->weight(1,0)*W;
234 
235  //Add the mass matrix term
236  if(flag==2)
237  {
238  mass_matrix(local_eqn,local_unknown)
239  += T[r]*test(l)*psi(l2)*W;
240  }
241 
242  //Add contribution to Elemental Matrix
243  for(unsigned i=0;i<DIM;i++)
244  {
245  //Temporary term used in assembly
246  double tmp = wind[i];
247  if(!ALE_is_disabled) tmp -= T[r]*mesh_velocity[i];
248  //Now assemble Jacobian term
249  jacobian(local_eqn,local_unknown)
250  -= dpsidx(l2,i)*(tmp*test(l)+D[r]*dtestdx(l,i))*W;
251  }
252 
253  } //End of diagonal terms
254 
255  //Now add the cross-reaction terms
256  jacobian(local_eqn,local_unknown) -=
257  dRdC(r,r2)*psi(l2)*test(l)*W;
258  }
259  }
260  }
261  } //End of jacobian
262  }
263  } //End of loop over reagents
264  } //End of loop over nodes
265  } // End of loop over integration points
266 }
267 
268 
269 
270 
271 //======================================================================
272 /// Self-test: Return 0 for OK
273 //======================================================================
274 template <unsigned NREAGENT, unsigned DIM>
276 {
277 
278  bool passed=true;
279 
280  // Check lower-level stuff
281  if (FiniteElement::self_test()!=0)
282  {
283  passed=false;
284  }
285 
286  // Return verdict
287  if (passed)
288  {
289  return 0;
290  }
291  else
292  {
293  return 1;
294  }
295 
296 }
297 
298 
299 //=========================================================================
300 /// Integrate the reagent concentrations over the element
301 //========================================================================
302 template <unsigned NREAGENT, unsigned DIM>
305 {
306  //Find out how many nodes there are
307  const unsigned n_node = nnode();
308 
309  //Get the nodal index at which the unknown is stored
310  unsigned c_nodal_index[NREAGENT];
311  for(unsigned r=0;r<NREAGENT;r++)
312  {c_nodal_index[r]= c_index_adv_diff_react(r);}
313 
314  //Set up memory for the shape and test functions
315  Shape psi(n_node);
316  DShape dpsidx(n_node,DIM);
317 
318  //Set the value of n_intpt
319  unsigned n_intpt = integral_pt()->nweight();
320 
321  //Loop over the integration points
322  for(unsigned ipt=0;ipt<n_intpt;ipt++)
323  {
324  //Get the integral weight
325  double w = integral_pt()->weight(ipt);
326 
327  //Call the derivatives of the shape and test functions
328  double J = dshape_eulerian_at_knot(ipt,psi,dpsidx);
329 
330  //Premultiply the weights and the Jacobian
331  double W = w*J;
332 
333  //Calculate local values of the solution and its derivatives
334  //Allocate
335  Vector<double> interpolated_c(NREAGENT,0.0);
336 
337  //Calculate function value and derivatives:
338  // Loop over nodes
339  for(unsigned l=0;l<n_node;l++)
340  {
341  //Loop over the unknown reagents
342  for(unsigned r=0;r<NREAGENT;r++)
343  {
344  //Get the value at the node
345  const double c_value = raw_nodal_value(l,c_nodal_index[r]);
346 
347  //Calculate the interpolated value
348  interpolated_c[r] += c_value*psi(l);
349  }
350  }
351 
352  for(unsigned r=0;r<NREAGENT;r++)
353  {
354  C[r] += interpolated_c[r]*W;
355  }
356  } //End of loop over integration points
357 
358 }
359 
360 
361 
362 
363 
364 
365 //======================================================================
366 /// \short Output function:
367 ///
368 /// x,y,u,w_x,w_y or x,y,z,u,w_x,w_y,w_z
369 ///
370 /// nplot points in each coordinate direction
371 //======================================================================
372 template <unsigned NREAGENT, unsigned DIM>
374 output(std::ostream &outfile, const unsigned &nplot)
375 {
376  //Vector of local coordinates
377  Vector<double> s(DIM);
378 
379 
380  // Tecplot header info
381  outfile << tecplot_zone_string(nplot);
382 
383  // Loop over plot points
384  unsigned num_plot_points=nplot_points(nplot);
385  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
386  {
387  // Get local coordinates of plot point
388  get_s_plot(iplot,nplot,s);
389 
390  // Get Eulerian coordinate of plot point
391  Vector<double> x(DIM);
392  interpolated_x(s,x);
393 
394  for(unsigned i=0;i<DIM;i++)
395  {
396  outfile << x[i] << " ";
397  }
398  for(unsigned i=0;i<NREAGENT;i++)
399  {
400  outfile << interpolated_c_adv_diff_react(s,i) << " ";
401  }
402 
403  // Get the wind
404  Vector<double> wind(DIM);
405  // Dummy integration point needed
406  unsigned ipt=0;
407  get_wind_adv_diff_react(ipt,s,x,wind);
408  for(unsigned i=0;i<DIM;i++)
409  {
410  outfile << wind[i] << " ";
411  }
412  outfile << std::endl;
413 
414  }
415 
416  // Write tecplot footer (e.g. FE connectivity lists)
417  write_tecplot_zone_footer(outfile,nplot);
418 
419 }
420 
421 
422 //======================================================================
423 /// C-style output function:
424 ///
425 /// x,y,u or x,y,z,u
426 ///
427 /// nplot points in each coordinate direction
428 //======================================================================
429 template <unsigned NREAGENT,unsigned DIM>
431 output(FILE* file_pt, const unsigned &nplot)
432 {
433  //Vector of local coordinates
434  Vector<double> s(DIM);
435 
436  // Tecplot header info
437  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
438 
439  // Loop over plot points
440  unsigned num_plot_points=nplot_points(nplot);
441  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
442  {
443 
444  // Get local coordinates of plot point
445  get_s_plot(iplot,nplot,s);
446 
447  for(unsigned i=0;i<DIM;i++)
448  {
449  fprintf(file_pt,"%g ",interpolated_x(s,i));
450 
451  }
452  for(unsigned i=0;i<NREAGENT;i++)
453  {
454  fprintf(file_pt,"%g \n",interpolated_c_adv_diff_react(s,i));
455  }
456  }
457 
458  // Write tecplot footer (e.g. FE connectivity lists)
459  write_tecplot_zone_footer(file_pt,nplot);
460 
461 }
462 
463 
464 
465 //======================================================================
466  /// \short Output exact solution
467  ///
468  /// Solution is provided via function pointer.
469  /// Plot at a given number of plot points.
470  ///
471  /// x,y,u_exact or x,y,z,u_exact
472 //======================================================================
473 template <unsigned NREAGENT,unsigned DIM>
475 output_fct(std::ostream &outfile,
476  const unsigned &nplot,
478 {
479 
480  //Vector of local coordinates
481  Vector<double> s(DIM);
482 
483  // Vector for coordintes
484  Vector<double> x(DIM);
485 
486  // Tecplot header info
487  outfile << tecplot_zone_string(nplot);
488 
489  // Exact solution Vector (here a scalar)
490  Vector<double> exact_soln(1);
491 
492  // Loop over plot points
493  unsigned num_plot_points=nplot_points(nplot);
494  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
495  {
496 
497  // Get local coordinates of plot point
498  get_s_plot(iplot,nplot,s);
499 
500  // Get x position as Vector
501  interpolated_x(s,x);
502 
503  // Get exact solution at this point
504  (*exact_soln_pt)(x,exact_soln);
505 
506  //Output x,y,...,u_exact
507  for(unsigned i=0;i<DIM;i++)
508  {
509  outfile << x[i] << " ";
510  }
511  outfile << exact_soln[0] << std::endl;
512  }
513 
514  // Write tecplot footer (e.g. FE connectivity lists)
515  write_tecplot_zone_footer(outfile,nplot);
516 
517 }
518 
519 
520 
521 
522 //======================================================================
523  /// \short Validate against exact solution
524  ///
525  /// Solution is provided via function pointer.
526  /// Plot error at a given number of plot points.
527  ///
528 //======================================================================
529 template <unsigned NREAGENT, unsigned DIM>
531 compute_error(std::ostream &outfile,
533  double& error, double& norm)
534 {
535 
536  // Initialise
537  error=0.0;
538  norm=0.0;
539 
540  //Vector of local coordinates
541  Vector<double> s(DIM);
542 
543  // Vector for coordintes
544  Vector<double> x(DIM);
545 
546  //Find out how many nodes there are in the element
547  unsigned n_node = nnode();
548 
549  Shape psi(n_node);
550 
551  //Set the value of n_intpt
552  unsigned n_intpt = integral_pt()->nweight();
553 
554  // Tecplot header info
555  outfile << "ZONE" << std::endl;
556 
557  // Exact solution Vector (here a scalar)
558  Vector<double> exact_soln(1);
559 
560  //Loop over the integration points
561  for(unsigned ipt=0;ipt<n_intpt;ipt++)
562  {
563 
564  //Assign values of s
565  for(unsigned i=0;i<DIM;i++)
566  {
567  s[i] = integral_pt()->knot(ipt,i);
568  }
569 
570  //Get the integral weight
571  double w = integral_pt()->weight(ipt);
572 
573  // Get jacobian of mapping
574  double J=J_eulerian(s);
575 
576  //Premultiply the weights and the Jacobian
577  double W = w*J;
578 
579  // Get x position as Vector
580  interpolated_x(s,x);
581 
582  // Get FE function value
583  double u_fe=interpolated_c_adv_diff_react(s,0);
584 
585  // Get exact solution at this point
586  (*exact_soln_pt)(x,exact_soln);
587 
588  //Output x,y,...,error
589  for(unsigned i=0;i<DIM;i++)
590  {
591  outfile << x[i] << " ";
592  }
593  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
594 
595  // Add to error and norm
596  norm+=exact_soln[0]*exact_soln[0]*W;
597  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
598 
599  }
600 }
601 
602 //====================================================================
603 // Force build of templates, only building binary reactions at present
604 //====================================================================
605 ///One reagent only (!)
609 
613 
617 
621 
622 //Two reagents
626 
630 
634 
638 
639 }
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
virtual void fill_in_generic_residual_contribution_adv_diff_react(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...
cstr elem_len * i
Definition: cfortran.h:607
QAdvectionDiffusionReactionElement elements are linear/quadrilateral/brick-shaped Advection Diffusion...
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
void output(std::ostream &outfile)
Output with default number of plot points.
void integrate_reagents(Vector< double > &C) const
Return the integrated reagent concentrations.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
A class for all elements that solve the Advection Diffusion Reaction equations using isoparametric el...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points.
static char t char * s
Definition: cfortran.h:572
static Vector< double > Default_dimensionless_number
Static default value for the dimensionless numbers.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720