MFEM  v3.0
ex8p.cpp
Go to the documentation of this file.
1 // MFEM Example 8 - Parallel Version
2 //
3 // Compile with: make ex8p
4 //
5 // Sample runs: mpirun -np 4 ex8p -m ../data/square-disc.mesh
6 // mpirun -np 4 ex8p -m ../data/star.mesh
7 // mpirun -np 4 ex8p -m ../data/escher.mesh
8 // mpirun -np 4 ex8p -m ../data/fichera.mesh
9 // mpirun -np 4 ex8p -m ../data/square-disc-p2.vtk
10 // mpirun -np 4 ex8p -m ../data/square-disc-p3.mesh
11 // mpirun -np 4 ex8p -m ../data/star-surf.mesh -o 2
12 //
13 // Description: This example code demonstrates the use of the Discontinuous
14 // Petrov-Galerkin (DPG) method in its primal 2x2 block form as a
15 // simple finite element discretization of the Laplace problem
16 // -Delta u = f with homogeneous Dirichlet boundary conditions. We
17 // use high-order continuous trial space, a high-order interfacial
18 // (trace) space, and a high-order discontinuous test space
19 // defining a local dual (H^{-1}) norm.
20 //
21 // The example highlights the use of interfacial (trace) finite
22 // elements and spaces, trace face integrators and the definition
23 // of block operators and preconditioners.
24 //
25 // We recommend viewing examples 1-5 before viewing this example.
26
27 #include "mfem.hpp"
28 #include <fstream>
29 #include <iostream>
30
31 using namespace std;
32 using namespace mfem;
33
34 int main(int argc, char *argv[])
35 {
36  // 1. Initialize MPI.
37  int num_procs, myid;
38  MPI_Init(&argc, &argv);
39  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
40  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
41
42  // 2. Parse command-line options.
43  const char *mesh_file = "../data/star.mesh";
44  int order = 1;
45  bool visualization = 1;
46
47  OptionsParser args(argc, argv);
49  "Mesh file to use.");
51  "Finite element order (polynomial degree).");
53  "--no-visualization",
54  "Enable or disable GLVis visualization.");
55  args.Parse();
56  if (!args.Good())
57  {
58  if (myid == 0)
59  args.PrintUsage(cout);
60  MPI_Finalize();
61  return 1;
62  }
63  if (myid == 0)
64  args.PrintOptions(cout);
65
66  // 3. Read the (serial) mesh from the given mesh file on all processors. We
67  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
68  // and volume meshes with the same code.
69  Mesh *mesh;
70  ifstream imesh(mesh_file);
71  if (!imesh)
72  {
73  if (myid == 0)
74  cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
75  MPI_Finalize();
76  return 2;
77  }
78  mesh = new Mesh(imesh, 1, 1);
79  imesh.close();
80  int dim = mesh->Dimension();
81
82  // 4. Refine the serial mesh on all processors to increase the resolution. In
83  // this example we do 'ref_levels' of uniform refinement. We choose
84  // 'ref_levels' to be the largest number that gives a final mesh with no
85  // more than 10,000 elements.
86  {
87  int ref_levels =
88  (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
89  for (int l = 0; l < ref_levels; l++)
90  mesh->UniformRefinement();
91  }
92
93  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
94  // this mesh further in parallel to increase the resolution. Once the
95  // parallel mesh is defined, the serial mesh can be deleted.
96  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
97  delete mesh;
98  {
99  int par_ref_levels = 1;
100  for (int l = 0; l < par_ref_levels; l++)
101  pmesh->UniformRefinement();
102  }
103
104  // 6. Define the trial, interfacial (trace) and test DPG spaces:
105  // - The trial space, x0_space, contains the non-interfacial unknowns and
106  // has the essential BC.
107  // - The interfacial space, xhat_space, contains the interfacial unknowns
108  // and does not have essential BC.
109  // - The test space, test_space, is an enriched space where the enrichment
110  // degree may depend on the spatial dimension of the domain.
111  int trial_order = order;
112  int trace_order = order - 1;
113  int test_order = order; // reduced order, full order is (order + dim - 1)
114  if (test_order < trial_order)
115  if(myid == 0)
116  cerr << "Warning, test space not enriched enough to handle primal"
117  << " trial space\n";
118
119  FiniteElementCollection *x0_fec, *xhat_fec, *test_fec;
120
121  x0_fec = new H1_FECollection(trial_order, dim);
122  xhat_fec = new RT_Trace_FECollection(trace_order, dim);
123  test_fec = new L2_FECollection(test_order, dim);
124
125  ParFiniteElementSpace *x0_space, *xhat_space, *test_space;
126
127  x0_space = new ParFiniteElementSpace(pmesh, x0_fec);
128  xhat_space = new ParFiniteElementSpace(pmesh, xhat_fec);
129  test_space = new ParFiniteElementSpace(pmesh, test_fec);
130
131  int glob_true_s0 = x0_space->GlobalTrueVSize();
132  int glob_true_s1 = xhat_space->GlobalTrueVSize();
133  int glob_true_s_test = test_space->GlobalTrueVSize();
134  if (myid == 0)
135  {
136  cout << "\nNumber of Unknowns:\n"
137  << " Trial space, X0 : " << glob_true_s0 << '\n'
138  << " Interface space, Xhat : " << glob_true_s1 << '\n'
139  << " Test space, Y : " << glob_true_s_test << "\n\n";
140  }
141
142  // 7. Set up the linear form F(.) which corresponds to the right-hand side of
143  // the FEM linear system, which in this case is (f,phi_i) where f=1.0 and
144  // phi_i are the basis functions in the test finite element fespace.
145  ConstantCoefficient one(1.0);
146  ParLinearForm * F = new ParLinearForm(test_space);
148  F->Assemble();
149
150  ParGridFunction * x0 = new ParGridFunction(x0_space);
151  *x0 = 0.;
152
153  // 8. Set up the mixed bilinear form for the primal trial unknowns, B0,
154  // the mixed bilinear form for the interfacial unknowns, Bhat,
155  // the inverse stiffness matrix on the discontinuous test space, Sinv,
156  // and the stiffness matrix on the continuous trial space, S0.
157  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
158  ess_bdr = 1;
159  Array<int> ess_dof;
160  x0_space->GetEssentialVDofs(ess_bdr, ess_dof);
161
162  ParMixedBilinearForm *B0 = new ParMixedBilinearForm(x0_space,test_space);
164  B0->Assemble();
165  B0->EliminateEssentialBCFromTrialDofs(ess_dof, *x0, *F);
166  B0->Finalize();
167
168  ParMixedBilinearForm *Bhat = new ParMixedBilinearForm(xhat_space,test_space);
170  Bhat->Assemble();
171  Bhat->Finalize();
172
173  ParBilinearForm *Sinv = new ParBilinearForm(test_space);
174  SumIntegrator *Sum = new SumIntegrator;
178  Sinv->Assemble();
179  Sinv->Finalize();
180
181  ParBilinearForm *S0 = new ParBilinearForm(x0_space);
184  S0->Assemble();
185  S0->EliminateEssentialBC(ess_bdr);
186  S0->Finalize();
187
188  HypreParMatrix * matB0 = B0->ParallelAssemble(); delete B0;
189  HypreParMatrix * matBhat = Bhat->ParallelAssemble(); delete Bhat;
190  HypreParMatrix * matSinv = Sinv->ParallelAssemble(); delete Sinv;
191  HypreParMatrix * matS0 = S0->ParallelAssemble(); delete S0;
192
193  // 9. Define the block structure of the problem, by creating the offset
194  // variables. Also allocate two BlockVector objects to store the solution
195  // and rhs.
196  enum {x0_var, xhat_var, NVAR};
197
198  int true_s0 = x0_space->TrueVSize();
199  int true_s1 = xhat_space->TrueVSize();
200  int true_s_test = test_space->TrueVSize();
201
202  Array<int> true_offsets(NVAR+1);
203  true_offsets[0] = 0;
204  true_offsets[1] = true_s0;
205  true_offsets[2] = true_s0+true_s1;
206
207  Array<int> true_offsets_test(2);
208  true_offsets_test[0] = 0;
209  true_offsets_test[1] = true_s_test;
210
211  BlockVector x(true_offsets), b(true_offsets);
212  x = 0.0;
213  b = 0.0;
214
215  // 10. Set up the 1x2 block Least Squares DPG operator, B = [B0 Bhat],
216  // the normal equation operator, A = B^t Sinv B, and
217  // the normal equation right-hand-size, b = B^t Sinv F.
218  BlockOperator B(true_offsets_test, true_offsets);
219  B.SetBlock(0, 0, matB0);
220  B.SetBlock(0, 1, matBhat);
221
222  RAPOperator A(B, *matSinv, B);
223
224  HypreParVector *trueF = F->ParallelAssemble();
225  {
226  HypreParVector SinvF(test_space);
227  matSinv->Mult(*trueF, SinvF);
228  B.MultTranspose(SinvF, b);
229  }
230
231  // 11. Set up a block-diagonal preconditioner for the 2x2 normal equation
232  //
233  // [ S0^{-1} 0 ]
234  // [ 0 Shat^{-1} ] Shat = (Bhat^T Sinv Bhat)
235  //
236  // corresponding to the primal (x0) and interfacial (xhat) unknowns.
237  HypreBoomerAMG *S0inv = new HypreBoomerAMG(*matS0);
238
239  HypreParMatrix *Shat = RAP(matSinv, matBhat);
240
241  HyprePCG *Shatinv = new HyprePCG(*Shat);
242  Shatinv->SetTol(1e-3);
243  Shatinv->SetMaxIter(200);
244  Shatinv->SetZeroInintialIterate();
245
246  BlockDiagonalPreconditioner P(true_offsets);
247  P.SetDiagonalBlock(0, S0inv);
248  P.SetDiagonalBlock(1, Shatinv);
249
250  // 12. Solve the normal equation system using the PCG iterative solver.
251  // Check the weighted norm of residual for the DPG least square problem.
252  // Wrap the primal variable in a GridFunction for visualization purposes.
253  CGSolver pcg(MPI_COMM_WORLD);
254  pcg.SetOperator(A);
255  pcg.SetPreconditioner(P);
256  pcg.SetRelTol(1e-6);
257  pcg.SetMaxIter(200);
258  pcg.SetPrintLevel(1);
259  pcg.Mult(b, x);
260
261  {
262  HypreParVector LSres(test_space), tmp(test_space);
263  B.Mult(x, LSres);
264  LSres -= *trueF;
265  matSinv->Mult(LSres, tmp);
266  double res = sqrt(InnerProduct(LSres, tmp));
267  if (myid == 0)
268  cout << "\n|| B0*x0 + Bhat*xhat - F ||_{S^-1} = " << res << endl;
269  }
270
271  x0->Distribute(x.GetBlock(x0_var));
272
273  // 13. Save the refined mesh and the solution in parallel. This output can
274  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
275  {
276  ostringstream mesh_name, sol_name;
277  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
278  sol_name << "sol." << setfill('0') << setw(6) << myid;
279
280  ofstream mesh_ofs(mesh_name.str().c_str());
281  mesh_ofs.precision(8);
282  pmesh->Print(mesh_ofs);
283
284  ofstream sol_ofs(sol_name.str().c_str());
285  sol_ofs.precision(8);
286  x0->Save(sol_ofs);
287  }
288
289  // 14. Send the solution by socket to a GLVis server.
290  if (visualization)
291  {
292  char vishost[] = "localhost";
293  int visport = 19916;
294  socketstream sol_sock(vishost, visport);
295  sol_sock << "parallel " << num_procs << " " << myid << "\n";
296  sol_sock.precision(8);
297  sol_sock << "solution\n" << *pmesh << *x0 << flush;
298  }
299
300  // 15. Free the used memory.
301  delete trueF;
302  delete Shatinv;
303  delete S0inv;
304  delete Shat;
305  delete matB0;
306  delete matBhat;
307  delete matSinv;
308  delete matS0;
309  delete x0;
310  delete F;
311  delete test_space;
312  delete xhat_space;
313  delete x0_space;
314  delete test_fec;
315  delete xhat_fec;
316  delete x0_fec;
317  delete pmesh;
318
319  MPI_Finalize();
320
321  return 0;
322 }
void SetTol(double tol)
Definition: hypre.cpp:1339
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:47
Definition: solvers.hpp:109
Integrator defining a sum of multiple Integrators.
Definition: bilininteg.hpp:139
Subclass constant coefficient.
Definition: coefficient.hpp:57
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:34
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:396
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:250
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Definition: bilininteg.hpp:149
void SetBlock(int iRow, int iCol, Operator *op)
Add a block op in the block-entry (iblock, jblock).
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: solvers.cpp:285
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:572
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void EliminateEssentialBCFromTrialDofs(Array< int > &marked_vdofs, Vector &sol, Vector &rhs)
Integrator that inverts the matrix assembled by another integrator.
Definition: bilininteg.hpp:121
The BoomerAMG solver in hypre.
Definition: hypre.hpp:519
Class for parallel linear form.
Definition: plinearform.hpp:26
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
void SetZeroInintialIterate()
non-hypre setting
Definition: hypre.hpp:404
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6225
void SetMaxIter(int max_it)
Definition: solvers.hpp:61
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:781
T Max() const
Definition: array.cpp:78
void Assemble(int skip_zeros=1)
void Assemble(int skip_zeros=1)
Assemble the local matrix.
The operator x -&gt; R*A*P*x.
Definition: operator.hpp:164
int Dimension() const
Definition: mesh.hpp:417
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:385
void SetMaxIter(int max_iter)
Definition: hypre.cpp:1344
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:38
Definition: linearform.cpp:19
Array< int > bdr_attributes
Definition: mesh.hpp:305
int main(int argc, char *argv[])
Definition: ex1.cpp:39
void SetRelTol(double rtol)
Definition: solvers.hpp:59
PCG solver in hypre.
Definition: hypre.hpp:381
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)
Definition: optparser.hpp:74
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:78
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:167
Class for parallel bilinear form.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:266
Class for parallel bilinear form.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:123
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:90
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:52
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs) const
Determine the boundary degrees of freedom.
Definition: pfespace.cpp:430
Class for parallel grid function.
Definition: pgridfunc.hpp:31
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:103
A class to handle Block systems in a matrix-free implementation.
Class for parallel meshes.
Definition: pmesh.hpp:27
void EliminateEssentialBC(Array< int > &bdr_attr_is_ess, Vector &sol, Vector &rhs, int d=0)
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:34
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:83
virtual void Print(std::ostream &out=std::cout) const
Definition: pmesh.cpp:2444
bool Good() const
Definition: optparser.hpp:120