MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex33p.cpp
Go to the documentation of this file.
1 // MFEM Example 33 - Parallel Version
2 //
3 // Compile with: make ex33p
4 //
5 // Sample runs: mpirun -np 4 ex33p -m ../data/square-disc.mesh -alpha 0.33 -o 2
6 // mpirun -np 4 ex33p -m ../data/square-disc.mesh -alpha 4.5 -o 3
7 // mpirun -np 4 ex33p -m ../data/star.mesh -alpha 1.4 -o 3
8 // mpirun -np 4 ex33p -m ../data/star.mesh -alpha 0.99 -o 3
9 // mpirun -np 4 ex33p -m ../data/inline-quad.mesh -alpha 0.5 -o 3
10 // mpirun -np 4 ex33p -m ../data/amr-quad.mesh -alpha 1.5 -o 3
11 // mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 -r 2
12 // mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4
13 // mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4
14 // mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5
15 //
16 // Verification runs:
17 // mpirun -np 4 ex33p -m ../data/inline-segment.mesh -ver -alpha 1.7 -o 2 -r 2
18 // mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 1.2 -o 2 -r 2
19 // mpirun -np 4 ex33p -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2
20 // mpirun -np 4 ex33p -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  // 0. Initialize MPI.
90  Mpi::Init(argc, argv);
91  int num_procs = Mpi::WorldSize();
92  int myid = Mpi::WorldRank();
93  Hypre::Init();
94 
95  // 1. Parse command-line options.
96  const char *mesh_file = "../data/star.mesh";
97  int order = 1;
98  int num_refs = 3;
99  double alpha = 0.5;
100  bool visualization = true;
101  bool verification = false;
102 
103  OptionsParser args(argc, argv);
104  args.AddOption(&mesh_file, "-m", "--mesh",
105  "Mesh file to use.");
106  args.AddOption(&order, "-o", "--order",
107  "Finite element order (polynomial degree) or -1 for"
108  " isoparametric space.");
109  args.AddOption(&num_refs, "-r", "--refs",
110  "Number of uniform refinements");
111  args.AddOption(&alpha, "-alpha", "--alpha",
112  "Fractional exponent");
113  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
114  "--no-visualization",
115  "Enable or disable GLVis visualization.");
116  args.AddOption(&verification, "-ver", "--verification", "-no-ver",
117  "--no-verification",
118  "Use sinusoidal function (f) for analytic comparison.");
119  args.Parse();
120  if (!args.Good())
121  {
122  args.PrintUsage(cout);
123  return 1;
124  }
125  if (Mpi::Root())
126  {
127  args.PrintOptions(cout);
128  }
129 
130  Array<double> coeffs, poles;
131  int progress_steps = 1;
132 
133  // 2. Compute the rational expansion coefficients that define the
134  // integer-order PDEs.
135  const int power_of_laplace = floor(alpha);
136  double exponent_to_approximate = alpha - power_of_laplace;
137  bool integer_order = false;
138  // Check if alpha is an integer or not.
139  if (abs(exponent_to_approximate) > 1e-12)
140  {
141  if (Mpi::Root())
142  {
143  mfem::out << "Approximating the fractional exponent "
144  << exponent_to_approximate
145  << endl;
146  }
147  ComputePartialFractionApproximation(exponent_to_approximate, coeffs,
148  poles);
149 
150  // If the example is build without LAPACK, the exponent_to_approximate
151  // might be modified by the function call above.
152  alpha = exponent_to_approximate + power_of_laplace;
153  }
154  else
155  {
156  integer_order = true;
157  if (Mpi::Root())
158  {
159  mfem::out << "Treating integer order PDE." << endl;
160  }
161  }
162 
163  // 3. Read the mesh from the given mesh file.
164  Mesh mesh(mesh_file, 1, 1);
165  int dim = mesh.Dimension();
166 
167  // 4. Refine the mesh to increase the resolution.
168  for (int i = 0; i < num_refs; i++)
169  {
170  mesh.UniformRefinement();
171  }
172  ParMesh pmesh(MPI_COMM_WORLD, mesh);
173  mesh.Clear();
174 
175  // 5. Define a finite element space on the mesh.
176  H1_FECollection fec(order, dim);
177  ParFiniteElementSpace fespace(&pmesh, &fec);
178  if (Mpi::Root())
179  {
180  cout << "Number of finite element unknowns: "
181  << fespace.GetTrueVSize() << endl;
182  }
183 
184  // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
185  Array<int> ess_tdof_list;
186  if (pmesh.bdr_attributes.Size())
187  {
188  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
189  ess_bdr = 1;
190  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
191  }
192 
193  // 7. Define diffusion coefficient, load, and solution GridFunction.
194  auto func = [&alpha](const Vector &x)
195  {
196  double val = 1.0;
197  for (int i=0; i<x.Size(); i++)
198  {
199  val *= sin(M_PI*x(i));
200  }
201  return pow(x.Size()*pow(M_PI,2), alpha) * val;
202  };
203  FunctionCoefficient f(func);
204  ConstantCoefficient one(1.0);
205  ParGridFunction u(&fespace);
206  ParGridFunction x(&fespace);
207  ParGridFunction g(&fespace);
208  u = 0.0;
209  x = 0.0;
210  g = 0.0;
211 
212  // 8. Prepare for visualization.
213  char vishost[] = "localhost";
214  int visport = 19916;
215 
216  // 9. Set up the linear form b(.) for integer-order PDE solves.
217  ParLinearForm b(&fespace);
218  if (verification)
219  {
220  // This statement is only relevant for the verification of the code. It
221  // uses a different f such that an analytic solution is known and easy
222  // to compare with the numerical one. The FPDE becomes:
223  // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2
224  // -> u(x,y) = sin(\pi x) sin(\pi y)
226  }
227  else
228  {
230  }
231  b.Assemble();
232 
233  // ------------------------------------------------------------------------
234  // 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f.
235  // ------------------------------------------------------------------------
236 
237  if (power_of_laplace > 0)
238  {
239  // 10.1 Compute Stiffnes Matrix
240  ParBilinearForm k(&fespace);
242  k.Assemble();
243 
244  // 10.2 Compute Mass Matrix
245  ParBilinearForm m(&fespace);
247  m.Assemble();
248  HypreParMatrix mass;
249  Array<int> empty;
250  m.FormSystemMatrix(empty, mass);
251 
252  // 10.3 Form the system of equations
253  Vector B, X;
254  OperatorPtr Op;
255  k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B);
256  HypreBoomerAMG prec;
257  prec.SetPrintLevel(-1);
258  CGSolver cg(MPI_COMM_WORLD);
259  cg.SetRelTol(1e-12);
260  cg.SetMaxIter(2000);
261  cg.SetPrintLevel(3);
262  cg.SetPreconditioner(prec);
263  cg.SetOperator(*Op);
264 
265  if (Mpi::Root())
266  {
267  mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace
268  << " ( f ) " << endl;
269  }
270  for (int i = 0; i < power_of_laplace; i++)
271  {
272  // 10.4 Solve the linear system Op X = B (N times).
273  cg.Mult(B, X);
274  // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step
275  if (i == power_of_laplace - 1)
276  {
277  // Needed for visualization and solution verification.
278  k.RecoverFEMSolution(X, b, g);
279  if (integer_order && verification)
280  {
281  // For an integer order PDE, g is also our solution u.
282  u+=g;
283  }
284  if (visualization)
285  {
286  socketstream fout;
287  ostringstream oss_f;
288  fout.open(vishost, visport);
289  fout.precision(8);
290  oss_f.str(""); oss_f.clear();
291  oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ "
292  << power_of_laplace
293  << " g = f";
294  fout << "parallel " << num_procs << " " << myid << "\n"
295  << "solution\n" << pmesh << g
296  << "window_title '" << oss_f.str() << "'" << flush;
297  }
298  }
299 
300  // 10.6 Prepare for next iteration (primal / dual space)
301  mass.Mult(X, B);
302  X.SetSubVectorComplement(ess_tdof_list,0.0);
303  }
304 
305  // 10.7 Extract solution for the next step. The b now corresponds to the
306  // function g in the PDE.
307  const SparseMatrix* rm = fespace.GetRestrictionMatrix();
308  rm->MultTranspose(B, b);
309  }
310 
311  // ------------------------------------------------------------------------
312  // 11. Solve the fractional PDE by solving M integer order PDEs and adding
313  // up the solutions.
314  // ------------------------------------------------------------------------
315  if (!integer_order)
316  {
317  // Setup visualization.
318  socketstream xout, uout;
319  ostringstream oss_x, oss_u;
320  if (visualization)
321  {
322  xout.open(vishost, visport);
323  xout.precision(8);
324  uout.open(vishost, visport);
325  uout.precision(8);
326  }
327  // Iterate over all expansion coefficient that contribute to the
328  // solution.
329  for (int i = 0; i < coeffs.Size(); i++)
330  {
331  if (Mpi::Root())
332  {
333  mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
334  << " u = " << coeffs[i] << " g " << endl;
335  }
336 
337  // 11.1 Reset GridFunction for integer-order PDE solve.
338  x = 0.0;
339 
340  // 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve.
341  ParBilinearForm a(&fespace);
343  ConstantCoefficient d_i(-poles[i]);
345  a.Assemble();
346 
347  // 11.3 Assemble the bilinear form and the corresponding linear system.
348  OperatorPtr A;
349  Vector B, X;
350  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
351 
352  // 11.4 Solve the linear system A X = B.
353  HypreBoomerAMG prec;
354  prec.SetPrintLevel(-1);
355 
356  CGSolver cg(MPI_COMM_WORLD);
357  cg.SetRelTol(1e-12);
358  cg.SetMaxIter(2000);
359  cg.SetPrintLevel(3);
360  cg.SetPreconditioner(prec);
361  cg.SetOperator(*A);
362  cg.Mult(B, X);
363 
364  // 11.5 Recover the solution as a finite element grid function.
365  a.RecoverFEMSolution(X, b, x);
366 
367  // 11.6 Accumulate integer-order PDE solutions.
368  x *= coeffs[i];
369  u += x;
370 
371  // 11.7 Send fractional PDE solution to a GLVis server.
372  if (visualization)
373  {
374  oss_x.str(""); oss_x.clear();
375  oss_x << "Step " << progress_steps
376  << ": Solution of PDE -Δ u + " << -poles[i]
377  << " u = " << coeffs[i] << " g";
378  xout << "parallel " << num_procs << " " << myid << "\n"
379  << "solution\n" << pmesh << x
380  << "window_title '" << oss_x.str() << "'" << flush;
381 
382  oss_u.str(""); oss_u.clear();
383  oss_u << "Step " << progress_steps + 1
384  << ": Solution of fractional PDE (-Δ)^" << alpha
385  << " u = f";
386  uout << "parallel " << num_procs << " " << myid << "\n"
387  << "solution\n" << pmesh << u
388  << "window_title '" << oss_u.str() << "'"
389  << flush;
390  }
391  }
392  }
393 
394  // ------------------------------------------------------------------------
395  // 12. (optional) Verify the solution.
396  // ------------------------------------------------------------------------
397  if (verification)
398  {
399  auto solution = [] (const Vector &x)
400  {
401  double val = 1.0;
402  for (int i=0; i<x.Size(); i++)
403  {
404  val *= sin(M_PI*x(i));
405  }
406  return val;
407  };
408  FunctionCoefficient sol(solution);
409  double l2_error = u.ComputeL2Error(sol);
410 
411  if (Mpi::Root())
412  {
413  string analytic_solution,expected_mesh;
414  switch (dim)
415  {
416  case 1:
417  analytic_solution = "sin(π x)";
418  expected_mesh = "inline_segment.mesh";
419  break;
420  case 2:
421  analytic_solution = "sin(π x) sin(π y)";
422  expected_mesh = "inline_quad.mesh";
423  break;
424  default:
425  analytic_solution = "sin(π x) sin(π y) sin(π z)";
426  expected_mesh = "inline_hex.mesh";
427  break;
428  }
429 
430  mfem::out << "\n" << string(80,'=')
431  << "\n\nSolution Verification in "<< dim << "D \n\n"
432  << "Analytic solution : " << analytic_solution << "\n"
433  << "Expected mesh : " << expected_mesh <<"\n"
434  << "Your mesh : " << mesh_file << "\n"
435  << "L2 error : " << l2_error << "\n\n"
436  << string(80,'=') << endl;
437  }
438  }
439 
440  return 0;
441 }
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
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1032
Conjugate gradient method.
Definition: solvers.hpp:465
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
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 Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:711
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:289
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1579
Class for parallel linear form.
Definition: plinearform.hpp:26
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
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[]
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1659
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: pgridfunc.hpp:281
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void Assemble(int skip_zeros=1)
Assemble the local matrix.
double analytic_solution(const Vector &x)
Definition: ex7.cpp:252
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
void AddDomainIntegrator(LinearFormIntegrator *lfi)
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
void SetRelTol(double rtol)
Definition: solvers.hpp:198
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 void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
double a
Definition: lissajous.cpp:41
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:389
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
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Class for parallel bilinear form.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:899
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1732
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
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
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150