MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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 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 // 1. Parse command-line options.
98 const char *mesh_file = "../data/star.mesh";
99 int order = 1;
100 int num_refs = 3;
101 real_t alpha = 0.5;
102 bool visualization = true;
103 bool verification = false;
104
105 OptionsParser args(argc, argv);
106 args.AddOption(&mesh_file, "-m", "--mesh",
107 "Mesh file to use.");
108 args.AddOption(&order, "-o", "--order",
109 "Finite element order (polynomial degree) or -1 for"
110 " isoparametric space.");
111 args.AddOption(&num_refs, "-r", "--refs",
112 "Number of uniform refinements");
113 args.AddOption(&alpha, "-alpha", "--alpha",
114 "Fractional exponent");
115 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
116 "--no-visualization",
117 "Enable or disable GLVis visualization.");
118 args.AddOption(&verification, "-ver", "--verification", "-no-ver",
119 "--no-verification",
120 "Use sinusoidal function (f) for manufactured "
121 "solution test.");
122 args.Parse();
123 if (!args.Good())
124 {
125 args.PrintUsage(cout);
126 return 1;
127 }
128 args.PrintOptions(cout);
129
130 Array<real_t> 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 = (int)floor(alpha);
136 real_t 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 mfem::out << "Approximating the fractional exponent "
142 << exponent_to_approximate
143 << endl;
144 ComputePartialFractionApproximation(exponent_to_approximate, coeffs,
145 poles);
146
147 // If the example is built without LAPACK, the exponent_to_approximate
148 // might be modified by the function call above.
149 alpha = exponent_to_approximate + power_of_laplace;
150 }
151 else
152 {
153 integer_order = true;
154 mfem::out << "Treating integer order PDE." << endl;
155 }
156
157 // 3. Read the mesh from the given mesh file.
158 Mesh mesh(mesh_file, 1, 1);
159 int dim = mesh.Dimension();
160
161 // 4. Refine the mesh to increase the resolution.
162 for (int i = 0; i < num_refs; i++)
163 {
164 mesh.UniformRefinement();
165 }
166
167 // 5. Define a finite element space on the mesh.
168 H1_FECollection fec(order, dim);
169 FiniteElementSpace fespace(&mesh, &fec);
170 cout << "Number of degrees of freedom: "
171 << fespace.GetTrueVSize() << endl;
172
173 // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
174 Array<int> ess_tdof_list;
175 if (mesh.bdr_attributes.Size())
176 {
177 Array<int> ess_bdr(mesh.bdr_attributes.Max());
178 ess_bdr = 1;
179 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
180 }
181
182 // 7. Define diffusion coefficient, load, and solution GridFunction.
183 auto func = [&alpha](const Vector &x)
184 {
185 real_t val = 1.0;
186 for (int i=0; i<x.Size(); i++)
187 {
188 val *= sin(M_PI*x(i));
189 }
190 return pow(x.Size()*pow(M_PI,2), alpha) * val;
191 };
193 ConstantCoefficient one(1.0);
194 GridFunction u(&fespace);
195 GridFunction x(&fespace);
196 GridFunction g(&fespace);
197 u = 0.0;
198 x = 0.0;
199 g = 0.0;
200
201 // 8. Prepare for visualization.
202 char vishost[] = "localhost";
203 int visport = 19916;
204
205 // 9. Set up the linear form b(.) for integer-order PDE solves.
206 LinearForm b(&fespace);
207 if (verification)
208 {
209 // This statement is only relevant for the verification of the code. It
210 // uses a different f such that an analytic solution is known and easy
211 // to compare with the numerical one. The FPDE becomes:
212 // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2
213 // -> u(x,y) = sin(\pi x) sin(\pi y)
214 b.AddDomainIntegrator(new DomainLFIntegrator(f));
215 }
216 else
217 {
218 b.AddDomainIntegrator(new DomainLFIntegrator(one));
219 }
220 b.Assemble();
221
222 // ------------------------------------------------------------------------
223 // 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f.
224 // ------------------------------------------------------------------------
225
226 if (power_of_laplace > 0)
227 {
228 // 10.1 Compute Stiffnes Matrix
229 BilinearForm k(&fespace);
231 k.Assemble();
232
233 // 10.2 Compute Mass Matrix
234 BilinearForm m(&fespace);
236 m.Assemble();
237 SparseMatrix mass;
238 Array<int> empty;
239 m.FormSystemMatrix(empty, mass);
240
241 // 10.3 Form the system of equations
242 Vector B, X;
243 OperatorPtr Op;
244 k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B);
245 GSSmoother M((SparseMatrix&)(*Op));
246
247 mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace
248 << " ( f ) " << endl;
249 for (int i = 0; i < power_of_laplace; i++)
250 {
251 // 10.4 Solve the linear system Op X = B (N times).
252 PCG(*Op, M, B, X, 3, 300, 1e-12, 0.0);
253
254 // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step
255 if (i == power_of_laplace - 1)
256 {
257 // Needed for visualization and solution verification.
258 k.RecoverFEMSolution(X, b, g);
259 if (integer_order && verification)
260 {
261 // For an integer order PDE, g is also our solution u.
262 u+=g;
263 }
264 if (visualization)
265 {
266 socketstream fout;
267 ostringstream oss_f;
268 fout.open(vishost, visport);
269 fout.precision(8);
270 oss_f.str(""); oss_f.clear();
271 oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ "
272 << power_of_laplace
273 << " g = f";
274 fout << "solution\n" << mesh << g
275 << "window_title '" << oss_f.str() << "'" << flush;
276 }
277 }
278
279 // 10.6 Prepare for next iteration (primal / dual space)
280 mass.Mult(X, B);
281 X.SetSubVectorComplement(ess_tdof_list,0.0);
282 }
283
284 // 10.7 Extract solution for the next step. The b now corresponds to the
285 // function g in the PDE.
286 const SparseMatrix * R = fespace.GetRestrictionMatrix();
287 if (R)
288 {
289 R->MultTranspose(B,b);
290 }
291 else
292 {
293 b = B;
294 }
295 }
296
297 // ------------------------------------------------------------------------
298 // 11. Solve the fractional PDE by solving M integer order PDEs and adding
299 // up the solutions.
300 // ------------------------------------------------------------------------
301 if (!integer_order)
302 {
303 // Setup visualization.
304 socketstream xout, uout;
305 ostringstream oss_x, oss_u;
306 if (visualization)
307 {
308 xout.open(vishost, visport);
309 xout.precision(8);
310 uout.open(vishost, visport);
311 uout.precision(8);
312 }
313 // Iterate over all expansion coefficient that contribute to the
314 // solution.
315 for (int i = 0; i < coeffs.Size(); i++)
316 {
317 mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
318 << " u = " << coeffs[i] << " g " << endl;
319
320
321 // 11.1 Reset GridFunction for integer-order PDE solve.
322 x = 0.0;
323
324 // 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve.
325 BilinearForm a(&fespace);
326 a.AddDomainIntegrator(new DiffusionIntegrator(one));
327 ConstantCoefficient d_i(-poles[i]);
328 a.AddDomainIntegrator(new MassIntegrator(d_i));
329 a.Assemble();
330
331 // 11.3 Assemble the bilinear form and the corresponding linear system.
332 OperatorPtr A;
333 Vector B, X;
334 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
335
336 // 11.4 Solve the linear system A X = B.
337 GSSmoother M((SparseMatrix&)(*A));
338
339 PCG(*A, M, B, X, 3, 300, 1e-12, 0.0);
340
341 // 11.5 Recover the solution as a finite element grid function.
342 a.RecoverFEMSolution(X, b, x);
343
344 // 11.6 Accumulate integer-order PDE solutions.
345 x *= coeffs[i];
346 u += x;
347
348 // 11.7 Send fractional PDE solution to a GLVis server.
349 if (visualization)
350 {
351 oss_x.str(""); oss_x.clear();
352 oss_x << "Step " << progress_steps
353 << ": Solution of PDE -Δ u + " << -poles[i]
354 << " u = " << coeffs[i] << " g";
355 xout << "solution\n" << mesh << x
356 << "window_title '" << oss_x.str() << "'" << flush;
357
358 oss_u.str(""); oss_u.clear();
359 oss_u << "Step " << progress_steps + 1
360 << ": Solution of fractional PDE (-Δ)^" << alpha
361 << " u = f";
362 uout << "solution\n" << mesh << u
363 << "window_title '" << oss_u.str() << "'"
364 << flush;
365 }
366 }
367 }
368
369 // ------------------------------------------------------------------------
370 // 12. (optional) Verify the solution.
371 // ------------------------------------------------------------------------
372 if (verification)
373 {
374 auto solution = [] (const Vector &x)
375 {
376 real_t val = 1.0;
377 for (int i=0; i<x.Size(); i++)
378 {
379 val *= sin(M_PI*x(i));
380 }
381 return val;
382 };
383 FunctionCoefficient sol(solution);
384 real_t l2_error = u.ComputeL2Error(sol);
385
386 string manufactured_solution,expected_mesh;
387 switch (dim)
388 {
389 case 1:
390 manufactured_solution = "sin(π x)";
391 expected_mesh = "inline_segment.mesh";
392 break;
393 case 2:
394 manufactured_solution = "sin(π x) sin(π y)";
395 expected_mesh = "inline_quad.mesh";
396 break;
397 default:
398 manufactured_solution = "sin(π x) sin(π y) sin(π z)";
399 expected_mesh = "inline_hex.mesh";
400 break;
401 }
402
403 mfem::out << "\n" << string(80,'=')
404 << "\n\nSolution Verification in "<< dim << "D \n\n"
405 << "Manufactured solution : " << manufactured_solution << "\n"
406 << "Expected mesh : " << expected_mesh <<"\n"
407 << "Your mesh : " << mesh_file << "\n"
408 << "L2 error : " << l2_error << "\n\n"
409 << string(80,'=') << endl;
410 }
411
412 return 0;
413}
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
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
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 FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
A coefficient that is constant across space and time.
Class for domain integration .
Definition lininteg.hpp:109
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.hpp:620
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:588
A general function coefficient.
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Vector with associated FE space and LinearFormIntegrators.
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
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
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.
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.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
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()
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
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:913
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[]