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