MFEM  v4.5.1 Finite element discretization library
ex33.cpp
Go to the documentation of this file.
1 // MFEM Example 33
2 //
3 // Compile with: make ex33
4 //
5 // Sample runs: ex33 -m ../data/square-disc.mesh -alpha 0.33 -o 2
6 // ex33 -m ../data/square-disc.mesh -alpha 4.5 -o 3
7 // ex33 -m ../data/star.mesh -alpha 1.4 -o 3
8 // ex33 -m ../data/star.mesh -alpha 0.99 -o 3
9 // ex33 -m ../data/inline-quad.mesh -alpha 0.5 -o 3
10 // ex33 -m ../data/amr-quad.mesh -alpha 1.5 -o 3
11 // ex33 -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3
12 // ex33 -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4
13 // ex33 -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4
14 // ex33 -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5
15 //
16 // Verification runs:
17 // ex33 -m ../data/inline-segment.mesh -ver -alpha 1.7 -o 2 -r 2
18 // ex33 -m ../data/inline-quad.mesh -ver -alpha 1.2 -o 2 -r 2
19 // ex33 -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2
20 // ex33 -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1
21 //
22 // Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i)
23 // for all alpha.
24 //
25 // Description:
26 //
27 // In this example we solve the following fractional PDE with MFEM:
28 //
29 // ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α,
30 //
31 // To solve this FPDE, we apply the operator ( - Δ )^(-N), where the integer
32 // N is given by floor(α). By doing so, we obtain
33 //
34 // ( - Δ )^(α-N) u = ( - Δ )^(-N) f in Ω, u = 0 on ∂Ω, 0 < α.
35 //
36 // We first compute the right hand side by solving the integer order PDE
37 //
38 // ( - Δ )^N g = f in Ω, g = ( - Δ )^k g = 0 on ∂Ω, k = 1,..,N-1
39 //
40 // The remaining FPDE is then given by
41 //
42 // ( - Δ )^(α-N) u = g in Ω, u = 0 on ∂Ω.
43 //
44 // We rely on a rational approximation [2] of the normal linear operator
45 // A^{-α + N}, where A = - Δ (with associated homogeneous boundary conditions)
46 // and (a-N) in (0,1). We approximate the operator
47 //
48 // A^{-α+N} ≈ Σ_{i=0}^M c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0,
49 //
50 // where I is the L2-identity operator and the coefficients c_i and d_i
51 // are generated offline to a prescribed accuracy in a pre-processing step.
52 // We use the triple-A algorithm [1] to generate the rational approximation
53 // that this partial fractional expansion derives from. We then solve M+1
54 // independent integer-order PDEs,
55 //
56 // A u_i + d_i u_i = c_i g in Ω, u_i = 0 on ∂Ω, i=0,...,M,
57 //
58 // using MFEM and sum u_i to arrive at an approximate solution of the FPDE
59 //
60 // u ≈ Σ_{i=0}^M u_i.
61 //
62 // (If alpha is an integer, we stop after the first PDE was solved.)
63 //
64 // References:
65 //
66 // [1] Nakatsukasa, Y., Sète, O., & Trefethen, L. N. (2018). The AAA algorithm
67 // for rational approximation. SIAM Journal on Scientific Computing, 40(3),
68 // A1494-A1522.
69 //
70 // [2] Harizanov, S., Lazarov, R., Margenov, S., Marinov, P., & Pasciak, J.
71 // (2020). Analysis of numerical methods for spectral fractional elliptic
72 // equations based on the best uniform rational approximation. Journal of
73 // Computational Physics, 408, 109285.
74 //
75
76 #include "mfem.hpp"
77 #include <fstream>
78 #include <iostream>
79 #include <math.h>
80 #include <string>
81
82 #include "ex33.hpp"
83
84 using namespace std;
85 using namespace mfem;
86
87 int main(int argc, char *argv[])
88 {
89  // 1. Parse command-line options.
90  const char *mesh_file = "../data/star.mesh";
91  int order = 1;
92  int num_refs = 3;
93  double alpha = 0.5;
94  bool visualization = true;
95  bool verification = false;
96
97  OptionsParser args(argc, argv);
98  args.AddOption(&mesh_file, "-m", "--mesh",
99  "Mesh file to use.");
100  args.AddOption(&order, "-o", "--order",
101  "Finite element order (polynomial degree) or -1 for"
102  " isoparametric space.");
103  args.AddOption(&num_refs, "-r", "--refs",
104  "Number of uniform refinements");
105  args.AddOption(&alpha, "-alpha", "--alpha",
106  "Fractional exponent");
107  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
108  "--no-visualization",
109  "Enable or disable GLVis visualization.");
110  args.AddOption(&verification, "-ver", "--verification", "-no-ver",
111  "--no-verification",
112  "Use sinusoidal function (f) for analytic comparison.");
113  args.Parse();
114  if (!args.Good())
115  {
116  args.PrintUsage(cout);
117  return 1;
118  }
119  args.PrintOptions(cout);
120
121  Array<double> coeffs, poles;
122  int progress_steps = 1;
123
124  // 2. Compute the rational expansion coefficients that define the
125  // integer-order PDEs.
126  const int power_of_laplace = floor(alpha);
127  double exponent_to_approximate = alpha - power_of_laplace;
128  bool integer_order = false;
129  // Check if alpha is an integer or not.
130  if (abs(exponent_to_approximate) > 1e-12)
131  {
132  mfem::out << "Approximating the fractional exponent "
133  << exponent_to_approximate
134  << endl;
135  ComputePartialFractionApproximation(exponent_to_approximate, coeffs,
136  poles);
137
138  // If the example is build without LAPACK, the exponent_to_approximate
139  // might be modified by the function call above.
140  alpha = exponent_to_approximate + power_of_laplace;
141  }
142  else
143  {
144  integer_order = true;
145  mfem::out << "Treating integer order PDE." << endl;
146  }
147
148  // 3. Read the mesh from the given mesh file.
149  Mesh mesh(mesh_file, 1, 1);
150  int dim = mesh.Dimension();
151
152  // 4. Refine the mesh to increase the resolution.
153  for (int i = 0; i < num_refs; i++)
154  {
155  mesh.UniformRefinement();
156  }
157
158  // 5. Define a finite element space on the mesh.
159  H1_FECollection fec(order, dim);
160  FiniteElementSpace fespace(&mesh, &fec);
161  cout << "Number of finite element unknowns: "
162  << fespace.GetTrueVSize() << endl;
163
164  // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
165  Array<int> ess_tdof_list;
166  if (mesh.bdr_attributes.Size())
167  {
168  Array<int> ess_bdr(mesh.bdr_attributes.Max());
169  ess_bdr = 1;
170  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
171  }
172
173  // 7. Define diffusion coefficient, load, and solution GridFunction.
174  auto func = [&alpha](const Vector &x)
175  {
176  double val = 1.0;
177  for (int i=0; i<x.Size(); i++)
178  {
179  val *= sin(M_PI*x(i));
180  }
181  return pow(x.Size()*pow(M_PI,2), alpha) * val;
182  };
183  FunctionCoefficient f(func);
184  ConstantCoefficient one(1.0);
185  GridFunction u(&fespace);
186  GridFunction x(&fespace);
187  GridFunction g(&fespace);
188  u = 0.0;
189  x = 0.0;
190  g = 0.0;
191
192  // 8. Prepare for visualization.
193  char vishost[] = "localhost";
194  int visport = 19916;
195
196  // 9. Set up the linear form b(.) for integer-order PDE solves.
197  LinearForm b(&fespace);
198  if (verification)
199  {
200  // This statement is only relevant for the verification of the code. It
201  // uses a different f such that an analytic solution is known and easy
202  // to compare with the numerical one. The FPDE becomes:
203  // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2
204  // -> u(x,y) = sin(\pi x) sin(\pi y)
206  }
207  else
208  {
210  }
211  b.Assemble();
212
213  // ------------------------------------------------------------------------
214  // 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f.
215  // ------------------------------------------------------------------------
216
217  if (power_of_laplace > 0)
218  {
219  // 10.1 Compute Stiffnes Matrix
220  BilinearForm k(&fespace);
222  k.Assemble();
223
224  // 10.2 Compute Mass Matrix
225  BilinearForm m(&fespace);
227  m.Assemble();
228  SparseMatrix mass;
229  Array<int> empty;
230  m.FormSystemMatrix(empty, mass);
231
232  // 10.3 Form the system of equations
233  Vector B, X;
234  OperatorPtr Op;
235  k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B);
236  GSSmoother M((SparseMatrix&)(*Op));
237
238  mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace
239  << " ( f ) " << endl;
240  for (int i = 0; i < power_of_laplace; i++)
241  {
242  // 10.4 Solve the linear system Op X = B (N times).
243  PCG(*Op, M, B, X, 3, 300, 1e-12, 0.0);
244
245  // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step
246  if (i == power_of_laplace - 1)
247  {
248  // Needed for visualization and solution verification.
249  k.RecoverFEMSolution(X, b, g);
250  if (integer_order && verification)
251  {
252  // For an integer order PDE, g is also our solution u.
253  u+=g;
254  }
255  if (visualization)
256  {
257  socketstream fout;
258  ostringstream oss_f;
259  fout.open(vishost, visport);
260  fout.precision(8);
261  oss_f.str(""); oss_f.clear();
262  oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ "
263  << power_of_laplace
264  << " g = f";
265  fout << "solution\n" << mesh << g
266  << "window_title '" << oss_f.str() << "'" << flush;
267  }
268  }
269
270  // 10.6 Prepare for next iteration (primal / dual space)
271  mass.Mult(X, B);
272  X.SetSubVectorComplement(ess_tdof_list,0.0);
273  }
274
275  // 10.7 Extract solution for the next step. The b now corresponds to the
276  // function g in the PDE.
277  const SparseMatrix * R = fespace.GetRestrictionMatrix();
278  if (R)
279  {
280  R->MultTranspose(B,b);
281  }
282  else
283  {
284  b = B;
285  }
286  }
287
288  // ------------------------------------------------------------------------
289  // 11. Solve the fractional PDE by solving M integer order PDEs and adding
290  // up the solutions.
291  // ------------------------------------------------------------------------
292  if (!integer_order)
293  {
294  // Setup visualization.
295  socketstream xout, uout;
296  ostringstream oss_x, oss_u;
297  if (visualization)
298  {
299  xout.open(vishost, visport);
300  xout.precision(8);
301  uout.open(vishost, visport);
302  uout.precision(8);
303  }
304  // Iterate over all expansion coefficient that contribute to the
305  // solution.
306  for (int i = 0; i < coeffs.Size(); i++)
307  {
308  mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
309  << " u = " << coeffs[i] << " g " << endl;
310
311
312  // 11.1 Reset GridFunction for integer-order PDE solve.
313  x = 0.0;
314
315  // 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve.
316  BilinearForm a(&fespace);
318  ConstantCoefficient d_i(-poles[i]);
320  a.Assemble();
321
322  // 11.3 Assemble the bilinear form and the corresponding linear system.
323  OperatorPtr A;
324  Vector B, X;
325  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
326
327  // 11.4 Solve the linear system A X = B.
328  GSSmoother M((SparseMatrix&)(*A));
329
330  PCG(*A, M, B, X, 3, 300, 1e-12, 0.0);
331
332  // 11.5 Recover the solution as a finite element grid function.
333  a.RecoverFEMSolution(X, b, x);
334
335  // 11.6 Accumulate integer-order PDE solutions.
336  x *= coeffs[i];
337  u += x;
338
339  // 11.7 Send fractional PDE solution to a GLVis server.
340  if (visualization)
341  {
342  oss_x.str(""); oss_x.clear();
343  oss_x << "Step " << progress_steps
344  << ": Solution of PDE -Δ u + " << -poles[i]
345  << " u = " << coeffs[i] << " g";
346  xout << "solution\n" << mesh << x
347  << "window_title '" << oss_x.str() << "'" << flush;
348
349  oss_u.str(""); oss_u.clear();
350  oss_u << "Step " << progress_steps + 1
351  << ": Solution of fractional PDE (-Δ)^" << alpha
352  << " u = f";
353  uout << "solution\n" << mesh << u
354  << "window_title '" << oss_u.str() << "'"
355  << flush;
356  }
357  }
358  }
359
360  // ------------------------------------------------------------------------
361  // 12. (optional) Verify the solution.
362  // ------------------------------------------------------------------------
363  if (verification)
364  {
365  auto solution = [] (const Vector &x)
366  {
367  double val = 1.0;
368  for (int i=0; i<x.Size(); i++)
369  {
370  val *= sin(M_PI*x(i));
371  }
372  return val;
373  };
374  FunctionCoefficient sol(solution);
375  double l2_error = u.ComputeL2Error(sol);
376
377  string analytic_solution,expected_mesh;
378  switch (dim)
379  {
380  case 1:
381  analytic_solution = "sin(π x)";
382  expected_mesh = "inline_segment.mesh";
383  break;
384  case 2:
385  analytic_solution = "sin(π x) sin(π y)";
386  expected_mesh = "inline_quad.mesh";
387  break;
388  default:
389  analytic_solution = "sin(π x) sin(π y) sin(π z)";
390  expected_mesh = "inline_hex.mesh";
391  break;
392  }
393
394  mfem::out << "\n" << string(80,'=')
395  << "\n\nSolution Verification in "<< dim << "D \n\n"
396  << "Analytic solution : " << analytic_solution << "\n"
397  << "Expected mesh : " << expected_mesh <<"\n"
398  << "Your mesh : " << mesh_file << "\n"
399  << "L2 error : " << l2_error << "\n\n"
400  << string(80,'=') << endl;
401  }
402
403  return 0;
404 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:168
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:564
Data type for Gauss-Seidel smoother of sparse matrix.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:50
constexpr char vishost[]
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
constexpr int visport
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
double analytic_solution(const Vector &x)
Definition: ex7.cpp:252
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:904
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the dofs Array to the given val.
Definition: vector.cpp:708
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:729
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
void ComputePartialFractionApproximation(double &alpha, Array< double > &coeffs, Array< double > &poles, double lmax=1000., double tol=1e-10, int npoints=1000, int max_order=100)
Definition: ex33.hpp:295
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: gridfunc.cpp:2766
double a
Definition: lissajous.cpp:41
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
int dim
Definition: ex24.cpp:53
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:904
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:495
int main()
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150