MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
extrapolator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_EXTRAPOLATOR_HPP
13 #define MFEM_EXTRAPOLATOR_HPP
14 
15 #include "mfem.hpp"
16 
17 namespace mfem
18 {
19 
20 class DiscreteUpwindLOSolver;
21 class FluxBasedFCT;
22 
24 {
25 private:
26  Array<bool> &active_zones;
27  ParBilinearForm &M, &K;
28  HypreParMatrix *K_mat;
29  const Vector &b;
30 
31  DiscreteUpwindLOSolver *lo_solver;
32  Vector *lumpedM;
33 
34  void ComputeElementsMinMax(const ParGridFunction &gf,
35  Vector &el_min, Vector &el_max) const;
36  void ComputeBounds(const ParFiniteElementSpace &pfes,
37  const Vector &el_min, const Vector &el_max,
38  Vector &dof_min, Vector &dof_max) const;
39  void ZeroOutInactiveZones(Vector &dx);
40 
41 public:
42  // HO is standard FE advection solve; LO is upwind diffusion.
44 
46  ParBilinearForm &Kbf, const Vector &rhs);
47 
49 
50  virtual void Mult(const Vector &x, Vector &dx) const;
51 };
52 
53 // Extrapolates through DG advection based on:
54 // [1] Aslam, "A Partial Differential Equation Approach to Multidimensional
55 // Extrapolation", JCP 193(1), 2004.
56 // [2] Bochkov, Gibou, "PDE-Based Multidimensional Extrapolation of Scalar
57 // Fields over Interfaces with Kinks and High Curvatures", SISC 42(4), 2020.
59 {
60 public:
63  int xtrap_degree = 1;
64  bool visualization = false;
65  int vis_steps = 5;
66 
68 
69  // The known values taken from elements where level_set > 0, and extrapolated
70  // to all other elements. The known values are not changed.
71  void Extrapolate(Coefficient &level_set, const ParGridFunction &input,
72  const double time_period, ParGridFunction &xtrap);
73 
74  // Errors in cut elements, given an exact solution.
75  void ComputeLocalErrors(Coefficient &level_set, const ParGridFunction &exact,
76  const ParGridFunction &xtrap,
77  double &err_L1, double &err_L2, double &err_LI);
78 
79 private:
80  void TimeLoop(ParGridFunction &sltn, ODESolver &ode_solver, double t_final,
81  double dt, int vis_x_pos, std::string vis_name);
82 };
83 
85 {
86 private:
87  const ParGridFunction &ls_gf;
88 
89 public:
91  VectorCoefficient(ls.ParFESpace()->GetMesh()->Dimension()), ls_gf(ls) { }
92 
94 
95  virtual void Eval(Vector &V, ElementTransformation &T,
96  const IntegrationPoint &ip)
97  {
98  Vector grad_ls(vdim), n(vdim);
99  ls_gf.GetGradient(T, grad_ls);
100  const double norm_grad = grad_ls.Norml2();
101  V = grad_ls;
102  if (norm_grad > 0.0) { V /= norm_grad; }
103 
104  // Since positive level set values correspond to the known region, we
105  // transport into the opposite direction of the gradient.
106  V *= -1;
107  }
108 };
109 
111 {
112 private:
113  const ParGridFunction &u_gf;
114  int comp;
115 
116 public:
117  GradComponentCoeff(const ParGridFunction &u, int c) : u_gf(u), comp(c) { }
118 
119  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
120  {
121  Vector grad_u(T.GetDimension());
122  u_gf.GetGradient(T, grad_u);
123  return grad_u(comp);
124  }
125 };
126 
128 {
129 private:
130  const ParGridFunction &u_gf;
131  VectorCoefficient &n_coeff;
132 
133 public:
135  : u_gf(u), n_coeff(n) { }
136 
137  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
138  {
139  const int dim = T.GetDimension();
140  Vector n(dim), grad_u(dim);
141  n_coeff.Eval(n, T, ip);
142  u_gf.GetGradient(T, grad_u);
143  return n * grad_u;
144  }
145 };
146 
148 {
149 private:
150  const ParGridFunction &du_dx, &du_dy;
151  VectorCoefficient &n_coeff;
152 
153 public:
155  const ParGridFunction &dy, VectorCoefficient &n)
156  : du_dx(dx), du_dy(dy), n_coeff(n) { }
157 
158  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
159  {
160  const int dim = T.GetDimension();
161  Vector n(dim), grad_u(dim);
162  n_coeff.Eval(n, T, ip);
163  grad_u(0) = du_dx.GetValue(T, ip);
164  grad_u(1) = du_dy.GetValue(T, ip);
165  return n * grad_u;
166  }
167 };
168 
170 {
171 public:
173  const Vector &Mlump);
174 
175  void CalcLOSolution(const Vector &u, const Vector &rhs, Vector &du) const;
176 
177  Array<int> &GetKmap() { return K_smap; }
178 
179 protected:
181  const SparseMatrix &K;
182  mutable SparseMatrix D;
183 
185  const Vector &M_lumped;
186 
187  void ComputeDiscreteUpwindMatrix() const;
189 };
190 
191 } // namespace mfem
192 
193 #endif
enum mfem::Extrapolator::XtrapType xtrap_type
GradComponentCoeff(const ParGridFunction &u, int c)
Base class for vector Coefficients that optionally depend on time and space.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Base abstract class for first order time dependent operators.
Definition: operator.hpp:289
int GetDimension() const
Return the topological dimension of the reference element.
Definition: eltrans.hpp:165
void ComputeDiscreteUpwindMatrix() const
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
LevelSetNormalGradCoeff(const ParGridFunction &ls)
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: pgridfunc.cpp:264
NormalGradCoeff(const ParGridFunction &u, VectorCoefficient &n)
Abstract parallel finite element space.
Definition: pfespace.hpp:28
DiscreteUpwindLOSolver(ParFiniteElementSpace &space, const SparseMatrix &adv, const Vector &Mlump)
NormalGradComponentCoeff(const ParGridFunction &dx, const ParGridFunction &dy, VectorCoefficient &n)
const SparseMatrix & K
void Extrapolate(Coefficient &level_set, const ParGridFunction &input, const double time_period, ParGridFunction &xtrap)
Data type sparse matrix.
Definition: sparsemat.hpp:50
AdvectionOper(Array< bool > &zones, ParBilinearForm &Mbf, ParBilinearForm &Kbf, const Vector &rhs)
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
void CalcLOSolution(const Vector &u, const Vector &rhs, Vector &du) const
void ComputeLocalErrors(Coefficient &level_set, const ParGridFunction &exact, const ParGridFunction &xtrap, double &err_L1, double &err_L2, double &err_LI)
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
virtual void Mult(const Vector &x, Vector &dx) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
string space
Class for integration point with weight.
Definition: intrules.hpp:25
Mesh * GetMesh(int type)
Definition: ex29.cpp:218
int dim
Definition: ex24.cpp:53
enum mfem::AdvectionOper::AdvectionMode adv_mode
Class for parallel bilinear form.
Vector data type.
Definition: vector.hpp:60
ParFiniteElementSpace & pfes
AdvectionOper::AdvectionMode advection_mode
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
void ApplyDiscreteUpwindMatrix(ParGridFunction &u, Vector &du) const
void GetGradient(ElementTransformation &tr, Vector &grad) const
Definition: gridfunc.cpp:1645