MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
extrapolate.cpp
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 // ----------------------------------------------
13 // Extrapolation Miniapp: PDE-based extrapolation
14 // ----------------------------------------------
15 //
16 // This miniapp extrapolates a finite element function from a set of elements
17 // (known values) to the rest of the domain. The set of elements that contains
18 // the known values is specified by the positive values of a level set
19 // Coefficient. The known values are not modified. The miniapp supports two
20 // PDE-based approaches [1, 2], both of which rely on solving a sequence of
21 // advection problems in the direction of the unknown parts of the domain.
22 // The extrapolation can be constant (1st order), linear (2nd order), or
23 // quadratic (3rd order). These formal orders hold for a limited band around
24 // the zero level set, see the given references for more info.
25 //
26 // [1] Aslam, "A Partial Differential Equation Approach to Multidimensional
27 // Extrapolation", JCP 193(1), 2004.
28 // [2] Bochkov, Gibou, "PDE-Based Multidimensional Extrapolation of Scalar
29 // Fields over Interfaces with Kinks and High Curvatures", SISC 42(4), 2020.
30 //
31 // Compile with: make extrapolate
32 //
33 // Sample runs:
34 // mpirun -np 4 extrapolate -m "../../data/inline-segment.mesh" -rs 6 -ed 2
35 // mpirun -np 4 extrapolate -rs 5 -p 0 -ed 2
36 // mpirun -np 4 extrapolate -rs 5 -p 1 -ed 2
37 // mpirun -np 4 extrapolate -rs 5 -p 1 -et 1 -ed 1 -dg 1
38 // mpirun -np 4 extrapolate -m "../../data/inline-hex.mesh" -ed 1 -rs 1
39 // mpirun -np 4 extrapolate -m "../../data/inline-hex.mesh" -p 1 -ed 1 -rs 1
40 
41 #include "extrapolator.hpp"
42 
43 using namespace std;
44 using namespace mfem;
45 
46 int problem = 0;
47 
48 double domainLS(const Vector &coord)
49 {
50  // Map from [0,1] to [-1,1].
51  const int dim = coord.Size();
52  const double x = coord(0)*2.0 - 1.0,
53  y = (dim > 1) ? coord(1)*2.0 - 1.0 : 0.0,
54  z = (dim > 2) ? coord(2)*2.0 - 1.0 : 0.0;
55 
56  switch (problem)
57  {
58  case 0:
59  {
60  // Sphere.
61  return 0.75 - sqrt(x*x + y*y + z*z + 1e-12);
62  }
63  case 1:
64  {
65  // Star.
66  MFEM_VERIFY(dim > 1, "Problem 1 is not applicable to 1D.");
67 
68  return 0.60 - sqrt(x*x + y*y + z*z + 1e-12) +
69  0.25 * (y*y*y*y*y + 5.0*x*x*x*x*y - 10.0*x*x*y*y*y) /
70  pow(x*x + y*y + z*z + 1e-12, 2.5) *
71  std::cos(0.5*M_PI * z / 0.6);
72  }
73  default: MFEM_ABORT("Bad option for --problem!"); return 0.0;
74  }
75 }
76 
77 double solution0(const Vector &coord)
78 {
79  // Map from [0,1] to [-1,1].
80  const int dim = coord.Size();
81  const double x = coord(0)*2.0 - 1.0 + 0.25,
82  y = (dim > 1) ? coord(1)*2.0 - 1.0 : 0.0,
83  z = (dim > 2) ? coord(2)*2.0 - 1.0 : 0.0;
84 
85  return std::cos(M_PI * x) * std::cos(M_PI * y) * std::cos(M_PI * z);
86 }
87 
88 void PrintNorm(int myid, Vector &v, std::string text)
89 {
90  double norm = v.Norml1();
91  MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
92  if (myid == 0)
93  {
94  std::cout << std::setprecision(12) << std::fixed
95  << text << norm << std::endl;
96  }
97 }
98 
99 void PrintIntegral(int myid, ParGridFunction &g, std::string text)
100 {
101  ConstantCoefficient zero(0.0);
102  double norm = g.ComputeL1Error(zero);
103  if (myid == 0)
104  {
105  std::cout << std::setprecision(12) << std::fixed
106  << text << norm << std::endl;
107  }
108 }
109 
110 int main(int argc, char *argv[])
111 {
112  // Initialize MPI and HYPRE.
113  Mpi::Init();
114  int myid = Mpi::WorldRank();;
115  Hypre::Init();
116 
117  // Parse command-line options.
118  const char *mesh_file = "../../data/inline-quad.mesh";
119  int rs_levels = 2;
120  Extrapolator::XtrapType ex_type = Extrapolator::ASLAM;
121  AdvectionOper::AdvectionMode dg_mode = AdvectionOper::HO;
122  int ex_degree = 1;
123  int order = 2;
124  double distance = 0.35;
125  bool vis_on = true;
126  int vis_steps_cnt = 50;
127 
128  OptionsParser args(argc, argv);
129  args.AddOption(&mesh_file, "-m", "--mesh",
130  "Mesh file to use.");
131  args.AddOption(&rs_levels, "-rs", "--refine-serial",
132  "Number of times to refine the mesh uniformly in serial.");
133  args.AddOption((int*)&ex_type, "-et", "--extrap-type",
134  "Extrapolation type: Aslam (0) or Bochkov (1).");
135  args.AddOption((int*)&dg_mode, "-dg", "--dg-mode",
136  "DG advection mode: 0 - Standard High-Order,\n\t"
137  " 1 - Low-Order Upwind Diffusion.");
138  args.AddOption(&ex_degree, "-ed", "--extrap-degree",
139  "Extrapolation degree: 0/1/2 for constant/linear/quadratic.");
140  args.AddOption(&order, "-o", "--order",
141  "Finite element order (polynomial degree) or -1 for"
142  " isoparametric space.");
143  args.AddOption(&distance, "-d", "--distance",
144  "Extrapolation distance.");
145  args.AddOption(&problem, "-p", "--problem",
146  "0 - 2D circle,\n\t"
147  "1 - 2D star");
148  args.AddOption(&vis_on, "-vis", "--visualization", "-no-vis",
149  "--no-visualization",
150  "Enable or disable GLVis visualization.");
151  args.AddOption(&vis_steps_cnt, "-vs", "--visualization-steps",
152  "Visualize every n-th timestep.");
153  args.Parse();
154  if (!args.Good())
155  {
156  if (myid == 0) { args.PrintUsage(cout); }
157  return 1;
158  }
159  if (myid == 0) { args.PrintOptions(cout); }
160 
161  // Refine the mesh and distribute.
162  Mesh mesh(mesh_file, 1, 1);
163  for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
164  ParMesh pmesh(MPI_COMM_WORLD, mesh);
165  mesh.Clear();
166  const int dim = pmesh.Dimension();
167 
168  // Input function.
169  L2_FECollection fec_L2(order, dim);
170  ParFiniteElementSpace pfes_L2(&pmesh, &fec_L2);
171  ParGridFunction u(&pfes_L2);
172  FunctionCoefficient u0_coeff(solution0);
173  u.ProjectCoefficient(u0_coeff);
174 
175  // Extrapolate.
176  Extrapolator xtrap;
177  xtrap.xtrap_type = ex_type;
178  xtrap.advection_mode = dg_mode;
179  xtrap.xtrap_degree = ex_degree;
180  xtrap.visualization = vis_on;
181  xtrap.vis_steps = vis_steps_cnt;
182  FunctionCoefficient ls_coeff(domainLS);
183  ParGridFunction ux(&pfes_L2);
184  xtrap.Extrapolate(ls_coeff, u, distance, ux);
185 
186  PrintNorm(myid, ux, "Solution l1 norm: ");
187  PrintIntegral(myid, ux, "Solution L1 norm: ");
188 
189  GridFunctionCoefficient u_exact_coeff(&u);
190  double err_L1 = ux.ComputeL1Error(u_exact_coeff),
191  err_L2 = ux.ComputeL2Error(u_exact_coeff);
192  if (myid == 0)
193  {
194  std::cout << "Global L1 error: " << err_L1 << std::endl
195  << "Global L2 error: " << err_L2 << std::endl;
196  }
197  double loc_error_L1, loc_error_L2, loc_error_LI;
198  xtrap.ComputeLocalErrors(ls_coeff, u, ux,
199  loc_error_L1, loc_error_L2, loc_error_LI);
200  if (myid == 0)
201  {
202  std::cout << "Local L1 error: " << loc_error_L1 << std::endl
203  << "Local L2 error: " << loc_error_L2 << std::endl
204  << "Local Li error: " << loc_error_LI << std::endl;
205  }
206 
207  // ParaView output.
208  ParGridFunction ls_gf(&pfes_L2);
209  ls_gf.ProjectCoefficient(ls_coeff);
210  ParaViewDataCollection dacol("ParaViewExtrapolate", &pmesh);
211  dacol.SetLevelsOfDetail(order);
212  dacol.RegisterField("Level Set Function", &ls_gf);
213  dacol.RegisterField("Extrapolated Solution", &ux);
214  dacol.SetTime(1.0);
215  dacol.SetCycle(1);
216  dacol.Save();
217 
218  return 0;
219 }
enum mfem::Extrapolator::XtrapType xtrap_type
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
Helper class for ParaView visualization data.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
void PrintIntegral(int myid, ParGridFunction &g, std::string text)
Definition: extrapolate.cpp:99
void Extrapolate(Coefficient &level_set, const ParGridFunction &input, const double time_period, ParGridFunction &xtrap)
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
double domainLS(const Vector &coord)
Definition: extrapolate.cpp:48
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
void ComputeLocalErrors(Coefficient &level_set, const ParGridFunction &exact, const ParGridFunction &xtrap, double &err_L1, double &err_L2, double &err_LI)
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetTime(double t)
Set physical time (for time-dependent simulations)
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
virtual double ComputeL1Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:266
double Norml1() const
Returns the l_1 norm of the vector.
Definition: vector.cpp:821
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
int problem
Definition: ex15.cpp:62
void PrintNorm(int myid, Vector &v, std::string text)
Definition: extrapolate.cpp:88
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:281
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
int dim
Definition: ex24.cpp:53
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:904
double solution0(const Vector &coord)
Definition: extrapolate.cpp:77
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
virtual void Save() override
AdvectionOper::AdvectionMode advection_mode
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:284
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150