MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
extrapolate.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
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
43using namespace std;
44using namespace mfem;
45
46int problem = 0;
47
48real_t domainLS(const Vector &coord)
49{
50 // Map from [0,1] to [-1,1].
51 const int dim = coord.Size();
52 const real_t 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
77real_t solution0(const Vector &coord)
78{
79 // Map from [0,1] to [-1,1].
80 const int dim = coord.Size();
81 const real_t 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
88void PrintNorm(int myid, Vector &v, std::string text)
89{
90 real_t norm = v.Norml1();
91 MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
92 MPI_COMM_WORLD);
93 if (myid == 0)
94 {
95 std::cout << std::setprecision(12) << std::fixed
96 << text << norm << std::endl;
97 }
98}
99
100void PrintIntegral(int myid, ParGridFunction &g, std::string text)
101{
102 ConstantCoefficient zero(0.0);
103 real_t norm = g.ComputeL1Error(zero);
104 if (myid == 0)
105 {
106 std::cout << std::setprecision(12) << std::fixed
107 << text << norm << std::endl;
108 }
109}
110
111int main(int argc, char *argv[])
112{
113 // Initialize MPI and HYPRE.
114 Mpi::Init();
115 int myid = Mpi::WorldRank();
116 Hypre::Init();
117
118 // Parse command-line options.
119 const char *mesh_file = "../../data/inline-quad.mesh";
120 int rs_levels = 2;
123 int ex_degree = 1;
124 int order = 2;
125 real_t distance = 0.35;
126 bool vis_on = true;
127 int vis_steps_cnt = 50;
128
129 OptionsParser args(argc, argv);
131 "Mesh file to use.");
133 "Number of times to refine the mesh uniformly in serial.");
135 "Extrapolation type: Aslam (0) or Bochkov (1).");
137 "DG advection mode: 0 - Standard High-Order,\n\t"
138 " 1 - Low-Order Upwind Diffusion.");
140 "Extrapolation degree: 0/1/2 for constant/linear/quadratic.");
142 "Finite element order (polynomial degree) or -1 for"
143 " isoparametric space.");
145 "Extrapolation distance.");
147 "0 - 2D circle,\n\t"
148 "1 - 2D star");
150 "--no-visualization",
151 "Enable or disable GLVis visualization.");
153 "Visualize every n-th timestep.");
154 args.Parse();
155 if (!args.Good())
156 {
157 if (myid == 0) { args.PrintUsage(cout); }
158 return 1;
159 }
160 if (myid == 0) { args.PrintOptions(cout); }
161
162 // Refine the mesh and distribute.
163 Mesh mesh(mesh_file, 1, 1);
164 for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
165 ParMesh pmesh(MPI_COMM_WORLD, mesh);
166 mesh.Clear();
167 const int dim = pmesh.Dimension();
168
169 // Input function.
170 L2_FECollection fec_L2(order, dim);
171 ParFiniteElementSpace pfes_L2(&pmesh, &fec_L2);
172 ParGridFunction u(&pfes_L2);
174 u.ProjectCoefficient(u0_coeff);
175
176 // Extrapolate.
177 Extrapolator xtrap;
178 xtrap.xtrap_type = ex_type;
180 xtrap.xtrap_degree = ex_degree;
181 xtrap.visualization = vis_on;
182 xtrap.vis_steps = vis_steps_cnt;
184 ParGridFunction ux(&pfes_L2);
185 xtrap.Extrapolate(ls_coeff, u, distance, ux);
186
187 PrintNorm(myid, ux, "Solution l1 norm: ");
188 PrintIntegral(myid, ux, "Solution L1 norm: ");
189
190 GridFunctionCoefficient u_exact_coeff(&u);
191 real_t err_L1 = ux.ComputeL1Error(u_exact_coeff),
192 err_L2 = ux.ComputeL2Error(u_exact_coeff);
193 if (myid == 0)
194 {
195 std::cout << "Global L1 error: " << err_L1 << std::endl
196 << "Global L2 error: " << err_L2 << std::endl;
197 }
198 real_t loc_error_L1, loc_error_L2, loc_error_LI;
199 xtrap.ComputeLocalErrors(ls_coeff, u, ux,
200 loc_error_L1, loc_error_L2, loc_error_LI);
201 if (myid == 0)
202 {
203 std::cout << "Local L1 error: " << loc_error_L1 << std::endl
204 << "Local L2 error: " << loc_error_L2 << std::endl
205 << "Local Li error: " << loc_error_LI << std::endl;
206 }
207
208 // ParaView output.
209 ParGridFunction ls_gf(&pfes_L2);
210 ls_gf.ProjectCoefficient(ls_coeff);
211 ParaViewDataCollection dacol("ParaViewExtrapolate", &pmesh);
212 dacol.SetLevelsOfDetail(order);
213 dacol.RegisterField("Level Set Function", &ls_gf);
214 dacol.RegisterField("Extrapolated Solution", &ux);
215 dacol.SetTime(1.0);
216 dacol.SetCycle(1);
217 dacol.Save();
218
219 return 0;
220}
A coefficient that is constant across space and time.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
void ComputeLocalErrors(Coefficient &level_set, const ParGridFunction &exact, const ParGridFunction &xtrap, real_t &err_L1, real_t &err_L2, real_t &err_LI)
enum mfem::Extrapolator::XtrapType xtrap_type
void Extrapolate(Coefficient &level_set, const ParGridFunction &input, const real_t time_period, ParGridFunction &xtrap)
A general function coefficient.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Mesh data type.
Definition mesh.hpp:56
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
Definition pfespace.hpp:29
Class for parallel grid function.
Definition pgridfunc.hpp:33
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
real_t ComputeL1Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const override
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for parallel meshes.
Definition pmesh.hpp:34
Helper class for ParaView visualization data.
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
Vector data type.
Definition vector.hpp:80
real_t Norml1() const
Returns the l_1 norm of the vector.
Definition vector.cpp:861
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
int dim
Definition ex24.cpp:53
void PrintNorm(int myid, Vector &v, std::string text)
real_t solution0(const Vector &coord)
int problem
void PrintIntegral(int myid, ParGridFunction &g, std::string text)
real_t domainLS(const Vector &coord)
int main()
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
Helper struct to convert a C++ type to an MPI type.