1 // MFEM Example 4 - Parallel Version
2 //
3 // Compile with: make ex4p
4 //
5 // Sample runs: mpirun -np 4 ex4p -m ../data/square-disc.mesh
6 // mpirun -np 4 ex4p -m ../data/star.mesh
7 // mpirun -np 4 ex4p -m ../data/beam-tet.mesh
8 // mpirun -np 4 ex4p -m ../data/beam-hex.mesh
9 // mpirun -np 4 ex4p -m ../data/beam-hex.mesh -o 2 -pa
10 // mpirun -np 4 ex4p -m ../data/escher.mesh -o 2 -sc
11 // mpirun -np 4 ex4p -m ../data/fichera.mesh -o 2 -hb
12 // mpirun -np 4 ex4p -m ../data/fichera-q2.vtk
13 // mpirun -np 4 ex4p -m ../data/fichera-q3.mesh -o 2 -sc
14 // mpirun -np 4 ex4p -m ../data/square-disc-nurbs.mesh -o 3
15 // mpirun -np 4 ex4p -m ../data/beam-hex-nurbs.mesh -o 3
16 // mpirun -np 4 ex4p -m ../data/periodic-square.mesh -no-bc
17 // mpirun -np 4 ex4p -m ../data/periodic-cube.mesh -no-bc
18 // mpirun -np 4 ex4p -m ../data/amr-quad.mesh
19 // mpirun -np 3 ex4p -m ../data/amr-quad.mesh -o 2 -hb
20 // mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -sc
21 // mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb
22 // mpirun -np 4 ex4p -m ../data/star-surf.mesh -o 3 -hb
23 //
24 // Device sample runs:
25 // mpirun -np 4 ex4p -m ../data/star.mesh -pa -d cuda
26 // mpirun -np 4 ex4p -m ../data/star.mesh -pa -d raja-cuda
27 // mpirun -np 4 ex4p -m ../data/star.mesh -pa -d raja-omp
28 // mpirun -np 4 ex4p -m ../data/beam-hex.mesh -pa -d cuda
29 //
30 // Description: This example code solves a simple 2D/3D H(div) diffusion
31 // problem corresponding to the second order definite equation
32 // -grad(alpha div F) + beta F = f with boundary condition F dot n
33 // = <given normal field>. Here, we use a given exact solution F
34 // and compute the corresponding r.h.s. f. We discretize with
35 // Raviart-Thomas finite elements.
36 //
37 // The example demonstrates the use of H(div) finite element
38 // spaces with the grad-div and H(div) vector finite element mass
39 // bilinear form, as well as the computation of discretization
40 // error when the exact solution is known. Bilinear form
41 // hybridization and static condensation are also illustrated.
42 //
43 // We recommend viewing examples 1-3 before viewing this example.
44
45 #include "mfem.hpp"
46 #include <fstream>
47 #include <iostream>
48
49 using namespace std;
50 using namespace mfem;
51
52 // Exact solution, F, and r.h.s., f. See below for implementation.
53 void F_exact(const Vector &, Vector &);
54 void f_exact(const Vector &, Vector &);
55 double freq = 1.0, kappa;
56
main(int argc,char * argv[])57 int main(int argc, char *argv[])
58 {
59 // 1. Initialize MPI.
60 int num_procs, myid;
61 MPI_Init(&argc, &argv);
62 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
63 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
64
65 // 2. Parse command-line options.
66 const char *mesh_file = "../data/star.mesh";
67 int order = 1;
68 bool set_bc = true;
69 bool static_cond = false;
70 bool hybridization = false;
71 bool pa = false;
72 const char *device_config = "cpu";
73 bool visualization = 1;
74
75 OptionsParser args(argc, argv);
76 args.AddOption(&mesh_file, "-m", "--mesh",
77 "Mesh file to use.");
78 args.AddOption(&order, "-o", "--order",
79 "Finite element order (polynomial degree).");
80 args.AddOption(&set_bc, "-bc", "--impose-bc", "-no-bc", "--dont-impose-bc",
81 "Impose or not essential boundary conditions.");
82 args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
83 " solution.");
84 args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
85 "--no-static-condensation", "Enable static condensation.");
86 args.AddOption(&hybridization, "-hb", "--hybridization", "-no-hb",
87 "--no-hybridization", "Enable hybridization.");
88 args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
89 "--no-partial-assembly", "Enable Partial Assembly.");
90 args.AddOption(&device_config, "-d", "--device",
91 "Device configuration string, see Device::Configure().");
92 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
93 "--no-visualization",
94 "Enable or disable GLVis visualization.");
95 args.Parse();
96 if (!args.Good())
97 {
98 if (myid == 0)
99 {
100 args.PrintUsage(cout);
101 }
102 MPI_Finalize();
103 return 1;
104 }
105 if (myid == 0)
106 {
107 args.PrintOptions(cout);
108 }
109 kappa = freq * M_PI;
110
111 // 3. Enable hardware devices such as GPUs, and programming models such as
112 // CUDA, OCCA, RAJA and OpenMP based on command line options.
113 Device device(device_config);
114 if (myid == 0) { device.Print(); }
115
116 // 4. Read the (serial) mesh from the given mesh file on all processors. We
117 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
118 // and volume, as well as periodic meshes with the same code.
119 Mesh *mesh = new Mesh(mesh_file, 1, 1);
120 int dim = mesh->Dimension();
121 int sdim = mesh->SpaceDimension();
122
123 // 5. Refine the serial mesh on all processors to increase the resolution. In
124 // this example we do 'ref_levels' of uniform refinement. We choose
125 // 'ref_levels' to be the largest number that gives a final mesh with no
126 // more than 1,000 elements.
127 {
128 int ref_levels =
129 (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
130 for (int l = 0; l < ref_levels; l++)
131 {
132 mesh->UniformRefinement();
133 }
134 }
135
136 // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
137 // this mesh further in parallel to increase the resolution. Once the
138 // parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
139 // meshes need to be reoriented before we can define high-order Nedelec
140 // spaces on them (this is needed in the ADS solver below).
141 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
142 delete mesh;
143 {
144 int par_ref_levels = 2;
145 for (int l = 0; l < par_ref_levels; l++)
146 {
147 pmesh->UniformRefinement();
148 }
149 }
150 pmesh->ReorientTetMesh();
151
152 // 7. Define a parallel finite element space on the parallel mesh. Here we
153 // use the Raviart-Thomas finite elements of the specified order.
154 FiniteElementCollection *fec = new RT_FECollection(order-1, dim);
155 ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
156 HYPRE_BigInt size = fespace->GlobalTrueVSize();
157 if (myid == 0)
158 {
159 cout << "Number of finite element unknowns: " << size << endl;
160 }
161
162 // 8. Determine the list of true (i.e. parallel conforming) essential
163 // boundary dofs. In this example, the boundary conditions are defined
164 // by marking all the boundary attributes from the mesh as essential
165 // (Dirichlet) and converting them to a list of true dofs.
166 Array<int> ess_tdof_list;
167 if (pmesh->bdr_attributes.Size())
168 {
169 Array<int> ess_bdr(pmesh->bdr_attributes.Max());
170 ess_bdr = set_bc ? 1 : 0;
171 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
172 }
173
174 // 9. Set up the parallel linear form b(.) which corresponds to the
175 // right-hand side of the FEM linear system, which in this case is
176 // (f,phi_i) where f is given by the function f_exact and phi_i are the
177 // basis functions in the finite element fespace.
178 VectorFunctionCoefficient f(sdim, f_exact);
179 ParLinearForm *b = new ParLinearForm(fespace);
180 b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
181 b->Assemble();
182
183 // 10. Define the solution vector x as a parallel finite element grid function
184 // corresponding to fespace. Initialize x by projecting the exact
185 // solution. Note that only values from the boundary faces will be used
186 // when eliminating the non-homogeneous boundary condition to modify the
187 // r.h.s. vector b.
188 ParGridFunction x(fespace);
189 VectorFunctionCoefficient F(sdim, F_exact);
190 x.ProjectCoefficient(F);
191
192 // 11. Set up the parallel bilinear form corresponding to the H(div)
193 // diffusion operator grad alpha div + beta I, by adding the div-div and
194 // the mass domain integrators.
195 Coefficient *alpha = new ConstantCoefficient(1.0);
196 Coefficient *beta = new ConstantCoefficient(1.0);
197 ParBilinearForm *a = new ParBilinearForm(fespace);
198 if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
199 a->AddDomainIntegrator(new DivDivIntegrator(*alpha));
200 a->AddDomainIntegrator(new VectorFEMassIntegrator(*beta));
201
202 // 12. Assemble the parallel bilinear form and the corresponding linear
203 // system, applying any necessary transformations such as: parallel
204 // assembly, eliminating boundary conditions, applying conforming
205 // constraints for non-conforming AMR, static condensation,
206 // hybridization, etc.
207 FiniteElementCollection *hfec = NULL;
208 ParFiniteElementSpace *hfes = NULL;
209 if (static_cond)
210 {
211 a->EnableStaticCondensation();
212 }
213 else if (hybridization)
214 {
215 hfec = new DG_Interface_FECollection(order-1, dim);
216 hfes = new ParFiniteElementSpace(pmesh, hfec);
217 a->EnableHybridization(hfes, new NormalTraceJumpIntegrator(),
218 ess_tdof_list);
219 }
220 a->Assemble();
221
222 OperatorPtr A;
223 Vector B, X;
224 a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
225
226 if (myid == 0 && !pa)
227 {
228 cout << "Size of linear system: "
229 << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
230 }
231
232 // 13. Define and apply a parallel PCG solver for A X = B with the 2D AMS or
233 // the 3D ADS preconditioners from hypre. If using hybridization, the
234 // system is preconditioned with hypre's BoomerAMG. In the partial
235 // assembly case, use Jacobi preconditioning.
236 Solver *prec = NULL;
237 CGSolver *pcg = new CGSolver(MPI_COMM_WORLD);
238 pcg->SetOperator(*A);
239 pcg->SetRelTol(1e-12);
240 pcg->SetMaxIter(2000);
241 pcg->SetPrintLevel(1);
242 if (hybridization) { prec = new HypreBoomerAMG(*A.As<HypreParMatrix>()); }
243 else if (pa) { prec = new OperatorJacobiSmoother(*a, ess_tdof_list); }
244 else
245 {
246 ParFiniteElementSpace *prec_fespace =
247 (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
248 if (dim == 2) { prec = new HypreAMS(*A.As<HypreParMatrix>(), prec_fespace); }
249 else { prec = new HypreADS(*A.As<HypreParMatrix>(), prec_fespace); }
250 }
251 pcg->SetPreconditioner(*prec);
252 pcg->Mult(B, X);
253
254 // 14. Recover the parallel grid function corresponding to X. This is the
255 // local finite element solution on each processor.
256 a->RecoverFEMSolution(X, *b, x);
257
258 // 15. Compute and print the L^2 norm of the error.
259 {
260 double err = x.ComputeL2Error(F);
261 if (myid == 0)
262 {
263 cout << "\n|| F_h - F ||_{L^2} = " << err << '\n' << endl;
264 }
265 }
266
267 // 16. Save the refined mesh and the solution in parallel. This output can
268 // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
269 {
270 ostringstream mesh_name, sol_name;
271 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
272 sol_name << "sol." << setfill('0') << setw(6) << myid;
273
274 ofstream mesh_ofs(mesh_name.str().c_str());
275 mesh_ofs.precision(8);
276 pmesh->Print(mesh_ofs);
277
278 ofstream sol_ofs(sol_name.str().c_str());
279 sol_ofs.precision(8);
280 x.Save(sol_ofs);
281 }
282
283 // 17. Send the solution by socket to a GLVis server.
284 if (visualization)
285 {
286 char vishost[] = "localhost";
287 int visport = 19916;
288 socketstream sol_sock(vishost, visport);
289 sol_sock << "parallel " << num_procs << " " << myid << "\n";
290 sol_sock.precision(8);
291 sol_sock << "solution\n" << *pmesh << x << flush;
292 }
293
294 // 18. Free the used memory.
295 delete pcg;
296 delete prec;
297 delete hfes;
298 delete hfec;
299 delete a;
300 delete alpha;
301 delete beta;
302 delete b;
303 delete fespace;
304 delete fec;
305 delete pmesh;
306
307 MPI_Finalize();
308
309 return 0;
310 }
311
312
313 // The exact solution (for non-surface meshes)
F_exact(const Vector & p,Vector & F)314 void F_exact(const Vector &p, Vector &F)
315 {
316 int dim = p.Size();
317
318 double x = p(0);
319 double y = p(1);
320 // double z = (dim == 3) ? p(2) : 0.0; // Uncomment if F is changed to depend on z
321
322 F(0) = cos(kappa*x)*sin(kappa*y);
323 F(1) = cos(kappa*y)*sin(kappa*x);
324 if (dim == 3)
325 {
326 F(2) = 0.0;
327 }
328 }
329
330 // The right hand side
f_exact(const Vector & p,Vector & f)331 void f_exact(const Vector &p, Vector &f)
332 {
333 int dim = p.Size();
334
335 double x = p(0);
336 double y = p(1);
337 // double z = (dim == 3) ? p(2) : 0.0; // Uncomment if f is changed to depend on z
338
339 double temp = 1 + 2*kappa*kappa;
340
341 f(0) = temp*cos(kappa*x)*sin(kappa*y);
342 f(1) = temp*cos(kappa*y)*sin(kappa*x);
343 if (dim == 3)
344 {
345 f(2) = 0;
346 }
347 }
348