MFEM  v4.5.1 Finite element discretization library
ex18.hpp
Go to the documentation of this file.
1 // MFEM Example 18 - Serial/Parallel Shared Code
2
3 #include "mfem.hpp"
4
5 using namespace std;
6 using namespace mfem;
7
8 // Problem definition
9 extern int problem;
10
11 // Maximum characteristic speed (updated by integrators)
12 extern double max_char_speed;
13
14 extern const int num_equation;
15 extern const double specific_heat_ratio;
16 extern const double gas_constant;
17
18 // Time-dependent operator for the right-hand side of the ODE representing the
19 // DG weak form.
21 {
22 private:
23  const int dim;
24
25  FiniteElementSpace &vfes;
26  Operator &A;
27  SparseMatrix &Aflux;
28  DenseTensor Me_inv;
29
30  mutable Vector state;
31  mutable DenseMatrix f;
32  mutable DenseTensor flux;
33  mutable Vector z;
34
35  void GetFlux(const DenseMatrix &state_, DenseTensor &flux_) const;
36
37 public:
39  Operator &A_, SparseMatrix &Aflux_);
40
41  virtual void Mult(const Vector &x, Vector &y) const;
42
43  virtual ~FE_Evolution() { }
44 };
45
46 // Implements a simple Rusanov flux
48 {
49 private:
50  Vector flux1;
51  Vector flux2;
52
53 public:
54  RiemannSolver();
55  double Eval(const Vector &state1, const Vector &state2,
56  const Vector &nor, Vector &flux);
57 };
58
59 // Interior face term: <F.n(u),[w]>
61 {
62 private:
63  RiemannSolver rsolver;
64  Vector shape1;
65  Vector shape2;
66  Vector funval1;
67  Vector funval2;
68  Vector nor;
69  Vector fluxN;
70
71 public:
72  FaceIntegrator(RiemannSolver &rsolver_, const int dim);
73
74  virtual void AssembleFaceVector(const FiniteElement &el1,
75  const FiniteElement &el2,
77  const Vector &elfun, Vector &elvect);
78 };
79
80 // Implementation of class FE_Evolution
82  Operator &A_, SparseMatrix &Aflux_)
83  : TimeDependentOperator(A_.Height()),
84  dim(vfes_.GetFE(0)->GetDim()),
85  vfes(vfes_),
86  A(A_),
87  Aflux(Aflux_),
88  Me_inv(vfes.GetFE(0)->GetDof(), vfes.GetFE(0)->GetDof(), vfes.GetNE()),
89  state(num_equation),
90  f(num_equation, dim),
91  flux(vfes.GetNDofs(), dim, num_equation),
92  z(A.Height())
93 {
94  // Standard local assembly and inversion for energy mass matrices.
95  const int dof = vfes.GetFE(0)->GetDof();
96  DenseMatrix Me(dof);
97  DenseMatrixInverse inv(&Me);
98  MassIntegrator mi;
99  for (int i = 0; i < vfes.GetNE(); i++)
100  {
101  mi.AssembleElementMatrix(*vfes.GetFE(i), *vfes.GetElementTransformation(i), Me);
102  inv.Factor();
103  inv.GetInverseMatrix(Me_inv(i));
104  }
105 }
106
107 void FE_Evolution::Mult(const Vector &x, Vector &y) const
108 {
109  // 0. Reset wavespeed computation before operator application.
110  max_char_speed = 0.;
111
112  // 1. Create the vector z with the face terms -<F.n(u), [w]>.
113  A.Mult(x, z);
114
115  // 2. Add the element terms.
116  // i. computing the flux approximately as a grid function by interpolating
117  // at the solution nodes.
118  // ii. multiplying this grid function by a (constant) mixed bilinear form for
119  // each of the num_equation, computing (F(u), grad(w)) for each equation.
120
121  DenseMatrix xmat(x.GetData(), vfes.GetNDofs(), num_equation);
122  GetFlux(xmat, flux);
123
124  for (int k = 0; k < num_equation; k++)
125  {
126  Vector fk(flux(k).GetData(), dim * vfes.GetNDofs());
127  Vector zk(z.GetData() + k * vfes.GetNDofs(), vfes.GetNDofs());
129  }
130
131  // 3. Multiply element-wise by the inverse mass matrices.
132  Vector zval;
133  Array<int> vdofs;
134  const int dof = vfes.GetFE(0)->GetDof();
135  DenseMatrix zmat, ymat(dof, num_equation);
136
137  for (int i = 0; i < vfes.GetNE(); i++)
138  {
139  // Return the vdofs ordered byNODES
140  vfes.GetElementVDofs(i, vdofs);
141  z.GetSubVector(vdofs, zval);
142  zmat.UseExternalData(zval.GetData(), dof, num_equation);
143  mfem::Mult(Me_inv(i), zmat, ymat);
144  y.SetSubVector(vdofs, ymat.GetData());
145  }
146 }
147
148 // Physicality check (at end)
149 bool StateIsPhysical(const Vector &state, const int dim);
150
151 // Pressure (EOS) computation
152 inline double ComputePressure(const Vector &state, int dim)
153 {
154  const double den = state(0);
155  const Vector den_vel(state.GetData() + 1, dim);
156  const double den_energy = state(1 + dim);
157
158  double den_vel2 = 0;
159  for (int d = 0; d < dim; d++) { den_vel2 += den_vel(d) * den_vel(d); }
160  den_vel2 /= den;
161
162  return (specific_heat_ratio - 1.0) * (den_energy - 0.5 * den_vel2);
163 }
164
165 // Compute the vector flux F(u)
166 void ComputeFlux(const Vector &state, int dim, DenseMatrix &flux)
167 {
168  const double den = state(0);
169  const Vector den_vel(state.GetData() + 1, dim);
170  const double den_energy = state(1 + dim);
171
172  MFEM_ASSERT(StateIsPhysical(state, dim), "");
173
174  const double pres = ComputePressure(state, dim);
175
176  for (int d = 0; d < dim; d++)
177  {
178  flux(0, d) = den_vel(d);
179  for (int i = 0; i < dim; i++)
180  {
181  flux(1+i, d) = den_vel(i) * den_vel(d) / den;
182  }
183  flux(1+d, d) += pres;
184  }
185
186  const double H = (den_energy + pres) / den;
187  for (int d = 0; d < dim; d++)
188  {
189  flux(1+dim, d) = den_vel(d) * H;
190  }
191 }
192
193 // Compute the scalar F(u).n
194 void ComputeFluxDotN(const Vector &state, const Vector &nor,
195  Vector &fluxN)
196 {
197  // NOTE: nor in general is not a unit normal
198  const int dim = nor.Size();
199  const double den = state(0);
200  const Vector den_vel(state.GetData() + 1, dim);
201  const double den_energy = state(1 + dim);
202
203  MFEM_ASSERT(StateIsPhysical(state, dim), "");
204
205  const double pres = ComputePressure(state, dim);
206
207  double den_velN = 0;
208  for (int d = 0; d < dim; d++) { den_velN += den_vel(d) * nor(d); }
209
210  fluxN(0) = den_velN;
211  for (int d = 0; d < dim; d++)
212  {
213  fluxN(1+d) = den_velN * den_vel(d) / den + pres * nor(d);
214  }
215
216  const double H = (den_energy + pres) / den;
217  fluxN(1 + dim) = den_velN * H;
218 }
219
220 // Compute the maximum characteristic speed.
221 inline double ComputeMaxCharSpeed(const Vector &state, const int dim)
222 {
223  const double den = state(0);
224  const Vector den_vel(state.GetData() + 1, dim);
225
226  double den_vel2 = 0;
227  for (int d = 0; d < dim; d++) { den_vel2 += den_vel(d) * den_vel(d); }
228  den_vel2 /= den;
229
230  const double pres = ComputePressure(state, dim);
231  const double sound = sqrt(specific_heat_ratio * pres / den);
232  const double vel = sqrt(den_vel2 / den);
233
234  return vel + sound;
235 }
236
237 // Compute the flux at solution nodes.
238 void FE_Evolution::GetFlux(const DenseMatrix &x_, DenseTensor &flux_) const
239 {
240  const int flux_dof = flux_.SizeI();
241  const int flux_dim = flux_.SizeJ();
242
243  for (int i = 0; i < flux_dof; i++)
244  {
245  for (int k = 0; k < num_equation; k++) { state(k) = x_(i, k); }
246  ComputeFlux(state, flux_dim, f);
247
248  for (int d = 0; d < flux_dim; d++)
249  {
250  for (int k = 0; k < num_equation; k++)
251  {
252  flux_(i, d, k) = f(k, d);
253  }
254  }
255
256  // Update max char speed
257  const double mcs = ComputeMaxCharSpeed(state, flux_dim);
258  if (mcs > max_char_speed) { max_char_speed = mcs; }
259  }
260 }
261
262 // Implementation of class RiemannSolver
264  flux1(num_equation),
265  flux2(num_equation) { }
266
267 double RiemannSolver::Eval(const Vector &state1, const Vector &state2,
268  const Vector &nor, Vector &flux)
269 {
270  // NOTE: nor in general is not a unit normal
271  const int dim = nor.Size();
272
273  MFEM_ASSERT(StateIsPhysical(state1, dim), "");
274  MFEM_ASSERT(StateIsPhysical(state2, dim), "");
275
276  const double maxE1 = ComputeMaxCharSpeed(state1, dim);
277  const double maxE2 = ComputeMaxCharSpeed(state2, dim);
278
279  const double maxE = max(maxE1, maxE2);
280
281  ComputeFluxDotN(state1, nor, flux1);
282  ComputeFluxDotN(state2, nor, flux2);
283
284  double normag = 0;
285  for (int i = 0; i < dim; i++)
286  {
287  normag += nor(i) * nor(i);
288  }
289  normag = sqrt(normag);
290
291  for (int i = 0; i < num_equation; i++)
292  {
293  flux(i) = 0.5 * (flux1(i) + flux2(i))
294  - 0.5 * maxE * (state2(i) - state1(i)) * normag;
295  }
296
297  return maxE;
298 }
299
300 // Implementation of class FaceIntegrator
302  rsolver(rsolver_),
303  funval1(num_equation),
304  funval2(num_equation),
305  nor(dim),
306  fluxN(num_equation) { }
307
309  const FiniteElement &el2,
311  const Vector &elfun, Vector &elvect)
312 {
313  // Compute the term <F.n(u),[w]> on the interior faces.
314  const int dof1 = el1.GetDof();
315  const int dof2 = el2.GetDof();
316
317  shape1.SetSize(dof1);
318  shape2.SetSize(dof2);
319
320  elvect.SetSize((dof1 + dof2) * num_equation);
321  elvect = 0.0;
322
323  DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation);
324  DenseMatrix elfun2_mat(elfun.GetData() + dof1 * num_equation, dof2,
325  num_equation);
326
327  DenseMatrix elvect1_mat(elvect.GetData(), dof1, num_equation);
328  DenseMatrix elvect2_mat(elvect.GetData() + dof1 * num_equation, dof2,
329  num_equation);
330
331  // Integration order calculation from DGTraceIntegrator
332  int intorder;
333  if (Tr.Elem2No >= 0)
334  intorder = (min(Tr.Elem1->OrderW(), Tr.Elem2->OrderW()) +
335  2*max(el1.GetOrder(), el2.GetOrder()));
336  else
337  {
338  intorder = Tr.Elem1->OrderW() + 2*el1.GetOrder();
339  }
340  if (el1.Space() == FunctionSpace::Pk)
341  {
342  intorder++;
343  }
344  const IntegrationRule *ir = &IntRules.Get(Tr.GetGeometryType(), intorder);
345
346  for (int i = 0; i < ir->GetNPoints(); i++)
347  {
348  const IntegrationPoint &ip = ir->IntPoint(i);
349
350  Tr.SetAllIntPoints(&ip); // set face and element int. points
351
352  // Calculate basis functions on both elements at the face
353  el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
354  el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
355
356  // Interpolate elfun at the point
357  elfun1_mat.MultTranspose(shape1, funval1);
358  elfun2_mat.MultTranspose(shape2, funval2);
359
360  // Get the normal vector and the flux on the face
361  CalcOrtho(Tr.Jacobian(), nor);
362  const double mcs = rsolver.Eval(funval1, funval2, nor, fluxN);
363
364  // Update max char speed
365  if (mcs > max_char_speed) { max_char_speed = mcs; }
366
367  fluxN *= ip.weight;
368  for (int k = 0; k < num_equation; k++)
369  {
370  for (int s = 0; s < dof1; s++)
371  {
372  elvect1_mat(s, k) -= fluxN(k) * shape1(s);
373  }
374  for (int s = 0; s < dof2; s++)
375  {
376  elvect2_mat(s, k) += fluxN(k) * shape2(s);
377  }
378  }
379  }
380 }
381
382 // Check that the state is physical - enabled in debug mode
383 bool StateIsPhysical(const Vector &state, const int dim)
384 {
385  const double den = state(0);
386  const Vector den_vel(state.GetData() + 1, dim);
387  const double den_energy = state(1 + dim);
388
389  if (den < 0)
390  {
391  cout << "Negative density: ";
392  for (int i = 0; i < state.Size(); i++)
393  {
394  cout << state(i) << " ";
395  }
396  cout << endl;
397  return false;
398  }
399  if (den_energy <= 0)
400  {
401  cout << "Negative energy: ";
402  for (int i = 0; i < state.Size(); i++)
403  {
404  cout << state(i) << " ";
405  }
406  cout << endl;
407  return false;
408  }
409
410  double den_vel2 = 0;
411  for (int i = 0; i < dim; i++) { den_vel2 += den_vel(i) * den_vel(i); }
412  den_vel2 /= den;
413
414  const double pres = (specific_heat_ratio - 1.0) * (den_energy - 0.5 * den_vel2);
415
416  if (pres <= 0)
417  {
418  cout << "Negative pressure: " << pres << ", state: ";
419  for (int i = 0; i < state.Size(); i++)
420  {
421  cout << state(i) << " ";
422  }
423  cout << endl;
424  return false;
425  }
426  return true;
427 }
428
429 // Initial condition
430 void InitialCondition(const Vector &x, Vector &y)
431 {
432  MFEM_ASSERT(x.Size() == 2, "");
433
434  double radius = 0, Minf = 0, beta = 0;
435  if (problem == 1)
436  {
437  // "Fast vortex"
439  Minf = 0.5;
440  beta = 1. / 5.;
441  }
442  else if (problem == 2)
443  {
444  // "Slow vortex"
446  Minf = 0.05;
447  beta = 1. / 50.;
448  }
449  else
450  {
451  mfem_error("Cannot recognize problem."
452  "Options are: 1 - fast vortex, 2 - slow vortex");
453  }
454
455  const double xc = 0.0, yc = 0.0;
456
457  // Nice units
458  const double vel_inf = 1.;
459  const double den_inf = 1.;
460
461  // Derive remainder of background state from this and Minf
462  const double pres_inf = (den_inf / specific_heat_ratio) * (vel_inf / Minf) *
463  (vel_inf / Minf);
464  const double temp_inf = pres_inf / (den_inf * gas_constant);
465
467  r2rad += (x(0) - xc) * (x(0) - xc);
468  r2rad += (x(1) - yc) * (x(1) - yc);
470
471  const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
472
473  const double velX = vel_inf * (1 - beta * (x(1) - yc) / radius * exp(
475  const double velY = vel_inf * beta * (x(0) - xc) / radius * exp(-0.5 * r2rad);
476  const double vel2 = velX * velX + velY * velY;
477
478  const double specific_heat = gas_constant * specific_heat_ratio * shrinv1;
479  const double temp = temp_inf - 0.5 * (vel_inf * beta) *
480  (vel_inf * beta) / specific_heat * exp(-r2rad);
481
482  const double den = den_inf * pow(temp/temp_inf, shrinv1);
483  const double pres = den * gas_constant * temp;
484  const double energy = shrinv1 * pres / den + 0.5 * vel2;
485
486  y(0) = den;
487  y(1) = den * velX;
488  y(2) = den * velY;
489  y(3) = den * energy;
490 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe_base.hpp:235
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:573
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:584
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:162
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
Base abstract class for first order time dependent operators.
Definition: operator.hpp:289
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:547
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void AddMult(const Vector &x, Vector &y, const double a=1.0) const
y += A * x (default) or y += a * A * x
Definition: sparsemat.cpp:736
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
FaceIntegrator(RiemannSolver &rsolver_, const int dim)
Definition: ex18.hpp:301
int Space() const
Returns the type of FunctionSpace on the element.
Definition: fe_base.hpp:337
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
const IntegrationPoint & GetElement2IntPoint()
Get a const reference to the integration point in neighboring element 2 corresponding to the currentl...
Definition: eltrans.hpp:587
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3890
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3878
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
ElementTransformation * Elem2
Definition: eltrans.hpp:522
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:209
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2654
double max_char_speed
Definition: ex18.cpp:59
RiemannSolver()
Definition: ex18.hpp:263
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
double Eval(const Vector &state1, const Vector &state2, const Vector &nor, Vector &flux)
Definition: ex18.hpp:267
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void ComputeFluxDotN(const Vector &state, const Vector &nor, Vector &fluxN)
Definition: ex18.hpp:194
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
Data type sparse matrix.
Definition: sparsemat.hpp:50
const int num_equation
Definition: ex18.cpp:54
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
void SetAllIntPoints(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition: eltrans.hpp:566
const double specific_heat_ratio
Definition: ex18.cpp:55
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator resulting from a face integral term...
Definition: ex18.hpp:308
void ComputeFlux(const Vector &state, int dim, DenseMatrix &flux)
Definition: ex18.hpp:166
int SizeI() const
Definition: densemat.hpp:999
double ComputeMaxCharSpeed(const Vector &state, const int dim)
Definition: ex18.hpp:221
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:647
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
int problem
Definition: ex15.cpp:62
const IntegrationPoint & GetElement1IntPoint()
Get a const reference to the integration point in neighboring element 1 corresponding to the currentl...
Definition: eltrans.hpp:577
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: ex18.hpp:107
int SizeJ() const
Definition: densemat.hpp:1000
FE_Evolution(FiniteElementSpace &vfes_, Operator &A_, SparseMatrix &Aflux_)
Definition: ex18.hpp:81
Class for integration point with weight.
Definition: intrules.hpp:25
Definition: distance.cpp:107
This class is used to express the local action of a general nonlinear finite element operator...
Definition: nonlininteg.hpp:27
int dim
Definition: ex24.cpp:53
double ComputePressure(const Vector &state, int dim)
Definition: ex18.hpp:152
void InitialCondition(const Vector &x, Vector &y)
Definition: ex18.hpp:430
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2781
ElementTransformation * Elem1
Definition: eltrans.hpp:522
virtual ~FE_Evolution()
Definition: ex18.hpp:43
Vector data type.
Definition: vector.hpp:60
RefCoord s[3]
bool StateIsPhysical(const Vector &state, const int dim)
Definition: ex18.hpp:383
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:80
Abstract operator.
Definition: operator.hpp:24
const double gas_constant
Definition: ex18.cpp:56
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:953