MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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 manufactured solution used in this problem is
23//
24// u = ∏_{i=0}^{dim-1} sin(π x_i) ,
25//
26// regardless of the value of alpha.
27//
28// Description:
29//
30// In this example we solve the following fractional PDE with MFEM:
31//
32// ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α,
33//
34// To solve this FPDE, we apply the operator ( - Δ )^(-N), where the integer
35// N is given by floor(α). By doing so, we obtain
36//
37// ( - Δ )^(α-N) u = ( - Δ )^(-N) f in Ω, u = 0 on ∂Ω, 0 < α.
38//
39// We first compute the right hand side by solving the integer order PDE
40//
41// ( - Δ )^N g = f in Ω, g = ( - Δ )^k g = 0 on ∂Ω, k = 1,..,N-1
42//
43// The remaining FPDE is then given by
44//
45// ( - Δ )^(α-N) u = g in Ω, u = 0 on ∂Ω.
46//
47// We rely on a rational approximation [2] of the normal linear operator
48// A^{-α + N}, where A = - Δ (with associated homogeneous boundary conditions)
49// and (a-N) in (0,1). We approximate the operator
50//
51// A^{-α+N} ≈ Σ_{i=0}^M c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0,
52//
53// where I is the L2-identity operator and the coefficients c_i and d_i
54// are generated offline to a prescribed accuracy in a pre-processing step.
55// We use the triple-A algorithm [1] to generate the rational approximation
56// that this partial fractional expansion derives from. We then solve M+1
57// independent integer-order PDEs,
58//
59// A u_i + d_i u_i = c_i g in Ω, u_i = 0 on ∂Ω, i=0,...,M,
60//
61// using MFEM and sum u_i to arrive at an approximate solution of the FPDE
62//
63// u ≈ Σ_{i=0}^M u_i.
64//
65// (If alpha is an integer, we stop after the first PDE was solved.)
66//
67// References:
68//
69// [1] Nakatsukasa, Y., Sète, O., & Trefethen, L. N. (2018). The AAA algorithm
70// for rational approximation. SIAM Journal on Scientific Computing, 40(3),
71// A1494-A1522.
72//
73// [2] Harizanov, S., Lazarov, R., Margenov, S., Marinov, P., & Pasciak, J.
74// (2020). Analysis of numerical methods for spectral fractional elliptic
75// equations based on the best uniform rational approximation. Journal of
76// Computational Physics, 408, 109285.
77//
78
79#include "mfem.hpp"
80#include <fstream>
81#include <iostream>
82#include <math.h>
83#include <string>
84
85#include "ex33.hpp"
86
87using namespace std;
88using namespace mfem;
89
90int main(int argc, char *argv[])
91{
92#ifdef MFEM_USE_SINGLE
93 cout << "This example is not supported in single precision.\n\n";
94 return MFEM_SKIP_RETURN_VALUE;
95#endif
96
97 // 0. Initialize MPI.
98 Mpi::Init(argc, argv);
99 int num_procs = Mpi::WorldSize();
100 int myid = Mpi::WorldRank();
101 Hypre::Init();
102
103 // 1. Parse command-line options.
104 const char *mesh_file = "../data/star.mesh";
105 int order = 1;
106 int num_refs = 3;
107 real_t alpha = 0.5;
108 bool visualization = true;
109 bool verification = false;
110
111 OptionsParser args(argc, argv);
112 args.AddOption(&mesh_file, "-m", "--mesh",
113 "Mesh file to use.");
114 args.AddOption(&order, "-o", "--order",
115 "Finite element order (polynomial degree) or -1 for"
116 " isoparametric space.");
117 args.AddOption(&num_refs, "-r", "--refs",
118 "Number of uniform refinements");
119 args.AddOption(&alpha, "-alpha", "--alpha",
120 "Fractional exponent");
121 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
122 "--no-visualization",
123 "Enable or disable GLVis visualization.");
124 args.AddOption(&verification, "-ver", "--verification", "-no-ver",
125 "--no-verification",
126 "Use sinusoidal function (f) for manufactured "
127 "solution test.");
128 args.Parse();
129 if (!args.Good())
130 {
131 args.PrintUsage(cout);
132 return 1;
133 }
134 if (Mpi::Root())
135 {
136 args.PrintOptions(cout);
137 }
138
139 Array<real_t> coeffs, poles;
140 int progress_steps = 1;
141
142 // 2. Compute the rational expansion coefficients that define the
143 // integer-order PDEs.
144 const int power_of_laplace = floor(alpha);
145 real_t exponent_to_approximate = alpha - power_of_laplace;
146 bool integer_order = false;
147 // Check if alpha is an integer or not.
148 if (abs(exponent_to_approximate) > 1e-12)
149 {
150 if (Mpi::Root())
151 {
152 mfem::out << "Approximating the fractional exponent "
153 << exponent_to_approximate
154 << endl;
155 }
156 ComputePartialFractionApproximation(exponent_to_approximate, coeffs,
157 poles);
158
159 // If the example is build without LAPACK, the exponent_to_approximate
160 // might be modified by the function call above.
161 alpha = exponent_to_approximate + power_of_laplace;
162 }
163 else
164 {
165 integer_order = true;
166 if (Mpi::Root())
167 {
168 mfem::out << "Treating integer order PDE." << endl;
169 }
170 }
171
172 // 3. Read the mesh from the given mesh file.
173 Mesh mesh(mesh_file, 1, 1);
174 int dim = mesh.Dimension();
175
176 // 4. Refine the mesh to increase the resolution.
177 for (int i = 0; i < num_refs; i++)
178 {
179 mesh.UniformRefinement();
180 }
181 ParMesh pmesh(MPI_COMM_WORLD, mesh);
182 mesh.Clear();
183
184 // 5. Define a finite element space on the mesh.
185 H1_FECollection fec(order, dim);
186 ParFiniteElementSpace fespace(&pmesh, &fec);
187 HYPRE_BigInt size = fespace.GlobalTrueVSize();
188 if (Mpi::Root())
189 {
190 cout << "Number of degrees of freedom: "
191 << size << endl;
192 }
193
194 // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
195 Array<int> ess_tdof_list;
196 if (pmesh.bdr_attributes.Size())
197 {
198 Array<int> ess_bdr(pmesh.bdr_attributes.Max());
199 ess_bdr = 1;
200 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
201 }
202
203 // 7. Define diffusion coefficient, load, and solution GridFunction.
204 auto func = [&alpha](const Vector &x)
205 {
206 real_t val = 1.0;
207 for (int i=0; i<x.Size(); i++)
208 {
209 val *= sin(M_PI*x(i));
210 }
211 return pow(x.Size()*pow(M_PI,2), alpha) * val;
212 };
214 ConstantCoefficient one(1.0);
215 ParGridFunction u(&fespace);
216 ParGridFunction x(&fespace);
217 ParGridFunction g(&fespace);
218 u = 0.0;
219 x = 0.0;
220 g = 0.0;
221
222 // 8. Prepare for visualization.
223 char vishost[] = "localhost";
224 int visport = 19916;
225
226 // 9. Set up the linear form b(.) for integer-order PDE solves.
227 ParLinearForm b(&fespace);
228 if (verification)
229 {
230 // This statement is only relevant for the verification of the code. It
231 // uses a different f such that an manufactured solution is known and easy
232 // to compare with the numerical one. The FPDE becomes:
233 // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2
234 // -> u(x,y) = sin(\pi x) sin(\pi y)
235 b.AddDomainIntegrator(new DomainLFIntegrator(f));
236 }
237 else
238 {
239 b.AddDomainIntegrator(new DomainLFIntegrator(one));
240 }
241 b.Assemble();
242
243 // ------------------------------------------------------------------------
244 // 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f.
245 // ------------------------------------------------------------------------
246
247 if (power_of_laplace > 0)
248 {
249 // 10.1 Compute Stiffnes Matrix
250 ParBilinearForm k(&fespace);
252 k.Assemble();
253
254 // 10.2 Compute Mass Matrix
255 ParBilinearForm m(&fespace);
257 m.Assemble();
258 HypreParMatrix mass;
259 Array<int> empty;
260 m.FormSystemMatrix(empty, mass);
261
262 // 10.3 Form the system of equations
263 Vector B, X;
264 OperatorPtr Op;
265 k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B);
266 HypreBoomerAMG prec;
267 prec.SetPrintLevel(-1);
268 CGSolver cg(MPI_COMM_WORLD);
269 cg.SetRelTol(1e-12);
270 cg.SetMaxIter(2000);
271 cg.SetPrintLevel(3);
272 cg.SetPreconditioner(prec);
273 cg.SetOperator(*Op);
274
275 if (Mpi::Root())
276 {
277 mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace
278 << " ( f ) " << endl;
279 }
280 for (int i = 0; i < power_of_laplace; i++)
281 {
282 // 10.4 Solve the linear system Op X = B (N times).
283 cg.Mult(B, X);
284 // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step
285 if (i == power_of_laplace - 1)
286 {
287 // Needed for visualization and solution verification.
288 k.RecoverFEMSolution(X, b, g);
289 if (integer_order && verification)
290 {
291 // For an integer order PDE, g is also our solution u.
292 u+=g;
293 }
294 if (visualization)
295 {
296 socketstream fout;
297 ostringstream oss_f;
298 fout.open(vishost, visport);
299 fout.precision(8);
300 oss_f.str(""); oss_f.clear();
301 oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ "
302 << power_of_laplace
303 << " g = f";
304 fout << "parallel " << num_procs << " " << myid << "\n"
305 << "solution\n" << pmesh << g
306 << "window_title '" << oss_f.str() << "'" << flush;
307 }
308 }
309
310 // 10.6 Prepare for next iteration (primal / dual space)
311 mass.Mult(X, B);
312 X.SetSubVectorComplement(ess_tdof_list,0.0);
313 }
314
315 // 10.7 Extract solution for the next step. The b now corresponds to the
316 // function g in the PDE.
317 const SparseMatrix* rm = fespace.GetRestrictionMatrix();
318 rm->MultTranspose(B, b);
319 }
320
321 // ------------------------------------------------------------------------
322 // 11. Solve the fractional PDE by solving M integer order PDEs and adding
323 // up the solutions.
324 // ------------------------------------------------------------------------
325 if (!integer_order)
326 {
327 // Setup visualization.
328 socketstream xout, uout;
329 ostringstream oss_x, oss_u;
330 if (visualization)
331 {
332 xout.open(vishost, visport);
333 xout.precision(8);
334 uout.open(vishost, visport);
335 uout.precision(8);
336 }
337 // Iterate over all expansion coefficient that contribute to the
338 // solution.
339 for (int i = 0; i < coeffs.Size(); i++)
340 {
341 if (Mpi::Root())
342 {
343 mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
344 << " u = " << coeffs[i] << " g " << endl;
345 }
346
347 // 11.1 Reset GridFunction for integer-order PDE solve.
348 x = 0.0;
349
350 // 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve.
351 ParBilinearForm a(&fespace);
352 a.AddDomainIntegrator(new DiffusionIntegrator(one));
353 ConstantCoefficient d_i(-poles[i]);
354 a.AddDomainIntegrator(new MassIntegrator(d_i));
355 a.Assemble();
356
357 // 11.3 Assemble the bilinear form and the corresponding linear system.
358 OperatorPtr A;
359 Vector B, X;
360 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
361
362 // 11.4 Solve the linear system A X = B.
363 HypreBoomerAMG prec;
364 prec.SetPrintLevel(-1);
365
366 CGSolver cg(MPI_COMM_WORLD);
367 cg.SetRelTol(1e-12);
368 cg.SetMaxIter(2000);
369 cg.SetPrintLevel(3);
370 cg.SetPreconditioner(prec);
371 cg.SetOperator(*A);
372 cg.Mult(B, X);
373
374 // 11.5 Recover the solution as a finite element grid function.
375 a.RecoverFEMSolution(X, b, x);
376
377 // 11.6 Accumulate integer-order PDE solutions.
378 x *= coeffs[i];
379 u += x;
380
381 // 11.7 Send fractional PDE solution to a GLVis server.
382 if (visualization)
383 {
384 oss_x.str(""); oss_x.clear();
385 oss_x << "Step " << progress_steps
386 << ": Solution of PDE -Δ u + " << -poles[i]
387 << " u = " << coeffs[i] << " g";
388 xout << "parallel " << num_procs << " " << myid << "\n"
389 << "solution\n" << pmesh << x
390 << "window_title '" << oss_x.str() << "'" << flush;
391
392 oss_u.str(""); oss_u.clear();
393 oss_u << "Step " << progress_steps + 1
394 << ": Solution of fractional PDE (-Δ)^" << alpha
395 << " u = f";
396 uout << "parallel " << num_procs << " " << myid << "\n"
397 << "solution\n" << pmesh << u
398 << "window_title '" << oss_u.str() << "'"
399 << flush;
400 }
401 }
402 }
403
404 // ------------------------------------------------------------------------
405 // 12. (optional) Verify the solution.
406 // ------------------------------------------------------------------------
407 if (verification)
408 {
409 auto solution = [] (const Vector &x)
410 {
411 real_t val = 1.0;
412 for (int i=0; i<x.Size(); i++)
413 {
414 val *= sin(M_PI*x(i));
415 }
416 return val;
417 };
418 FunctionCoefficient sol(solution);
419 real_t l2_error = u.ComputeL2Error(sol);
420
421 if (Mpi::Root())
422 {
423 string manufactured_solution,expected_mesh;
424 switch (dim)
425 {
426 case 1:
427 manufactured_solution = "sin(π x)";
428 expected_mesh = "inline_segment.mesh";
429 break;
430 case 2:
431 manufactured_solution = "sin(π x) sin(π y)";
432 expected_mesh = "inline_quad.mesh";
433 break;
434 default:
435 manufactured_solution = "sin(π x) sin(π y) sin(π z)";
436 expected_mesh = "inline_hex.mesh";
437 break;
438 }
439
440 mfem::out << "\n" << string(80,'=')
441 << "\n\nSolution Verification in "<< dim << "D \n\n"
442 << "Manufactured solution : " << manufactured_solution << "\n"
443 << "Expected mesh : " << expected_mesh <<"\n"
444 << "Your mesh : " << mesh_file << "\n"
445 << "L2 error : " << l2_error << "\n\n"
446 << string(80,'=') << endl;
447 }
448 }
449
450 return 0;
451}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Conjugate gradient method.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
A coefficient that is constant across space and time.
Class for domain integration .
Definition lininteg.hpp:109
A general function coefficient.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, real_t alpha=1.0, real_t beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition hypre.cpp:1815
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
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 bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of 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).
Pointer to an Operator of a specified type.
Definition handle.hpp:34
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.
Class for parallel bilinear form.
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(....
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
const SparseMatrix * GetRestrictionMatrix() const override
Get the R matrix which restricts a local dof vector to true dof vector.
Definition pfespace.hpp:389
Class for parallel grid function.
Definition pgridfunc.hpp:33
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
Data type sparse matrix.
Definition sparsemat.hpp:51
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSubVectorComplement(const Array< int > &dofs, const real_t val)
Set all vector entries NOT in the dofs Array to the given val.
Definition vector.cpp:739
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
void ComputePartialFractionApproximation(real_t &alpha, Array< real_t > &coeffs, Array< real_t > &poles, real_t lmax=1000., real_t tol=1e-10, int npoints=1000, int max_order=100)
Definition ex33.hpp:295
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
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
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]