MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex1p.cpp
Go to the documentation of this file.
1 // MFEM Example 1 - Parallel Version
2 // Caliper Modification
3 //
4 // Compile with: make ex1p
5 //
6 // Sample runs: mpirun -np 4 ex1p -m ../data/square-disc.mesh
7 // mpirun -np 4 ex1p -m ../data/star.mesh
8 // mpirun -np 4 ex1p -m ../data/star-mixed.mesh
9 // mpirun -np 4 ex1p -m ../data/escher.mesh
10 // mpirun -np 4 ex1p -m ../data/fichera.mesh
11 // mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
12 // mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
13 // mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
14 // mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
15 // mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
16 // mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
17 // mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
18 // mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
19 // mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
20 // mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
21 // mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
22 // mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
23 // mpirun -np 4 ex1p -m ../data/star-surf.mesh
24 // mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
25 // mpirun -np 4 ex1p -m ../data/inline-segment.mesh
26 // mpirun -np 4 ex1p -m ../data/amr-quad.mesh
27 // mpirun -np 4 ex1p -m ../data/amr-hex.mesh
28 // mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
29 // mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
30 //
31 // Device sample runs:
32 // mpirun -np 4 ex1p -pa -d cuda
33 // mpirun -np 4 ex1p -pa -d occa-cuda
34 // mpirun -np 4 ex1p -pa -d raja-omp
35 // mpirun -np 4 ex1p -pa -d ceed-cpu
36 // * mpirun -np 4 ex1p -pa -d ceed-cuda
37 // mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
38 // mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
39 //
40 // Description: This example is a copy of Example 1 instrumented with the
41 // Caliper performance profilinh library. Any option supported by
42 // the Caliper ConfigManager can be passed to the code using a
43 // configuration string after -p or --caliper flag. For more
44 // information, see the Caliper documentation.
45 //
46 // Examples: mpirun -np 4 ex1p --caliper runtime-report
47 // mpirun -np 4 ex1p --caliper runtime-report,mem.highwatermark,mpi-report
48 //
49 // The first run will return the default report. The second run will also output
50 // the memory high-water mark and time spent in MPI routines.
51 
52 #include "mfem.hpp"
53 #include <fstream>
54 #include <iostream>
55 
56 using namespace std;
57 using namespace mfem;
58 
59 int main(int argc, char *argv[])
60 {
61  // 1. Initialize MPI and HYPRE.
62  Mpi::Init(argc, argv);
63  int num_procs = Mpi::WorldSize();
64  int myid = Mpi::WorldRank();
65  Hypre::Init();
66  // Define Caliper ConfigManager
67  cali::ConfigManager mgr;
68  // Caliper instrumentation
69  MFEM_PERF_FUNCTION;
70 
71  // 2. Parse command-line options.
72  const char *mesh_file = "../../data/star.mesh";
73  int order = 1;
74  bool static_cond = false;
75  bool pa = false;
76  const char *device_config = "cpu";
77  bool visualization = true;
78  const char* cali_config = "runtime-report";
79 
80  OptionsParser args(argc, argv);
81  args.AddOption(&mesh_file, "-m", "--mesh",
82  "Mesh file to use.");
83  args.AddOption(&order, "-o", "--order",
84  "Finite element order (polynomial degree) or -1 for"
85  " isoparametric space.");
86  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
87  "--no-static-condensation", "Enable static condensation.");
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.AddOption(&cali_config, "-p", "--caliper",
96  "Caliper configuration string.");
97  args.Parse();
98  if (!args.Good())
99  {
100  if (myid == 0)
101  {
102  args.PrintUsage(cout);
103  }
104  return 1;
105  }
106  if (myid == 0)
107  {
108  args.PrintOptions(cout);
109  }
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  // Caliper configuration
117  mgr.add(cali_config);
118  mgr.start();
119 
120  // 4. Read the (serial) mesh from the given mesh file on all processors. We
121  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
122  // and volume meshes with the same code.
123  Mesh mesh(mesh_file, 1, 1);
124  int dim = mesh.Dimension();
125 
126  // 5. Refine the serial mesh on all processors to increase the resolution. In
127  // this example we do 'ref_levels' of uniform refinement. We choose
128  // 'ref_levels' to be the largest number that gives a final mesh with no
129  // more than 10,000 elements.
130  {
131  int ref_levels =
132  (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
133  for (int l = 0; l < ref_levels; l++)
134  {
135  mesh.UniformRefinement();
136  }
137  }
138 
139  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
140  // this mesh further in parallel to increase the resolution. Once the
141  // parallel mesh is defined, the serial mesh can be deleted.
142  ParMesh pmesh(MPI_COMM_WORLD, mesh);
143  mesh.Clear();
144  {
145  int par_ref_levels = 2;
146  for (int l = 0; l < par_ref_levels; l++)
147  {
148  pmesh.UniformRefinement();
149  }
150  }
151 
152  // 7. Define a parallel finite element space on the parallel mesh. Here we
153  // use continuous Lagrange finite elements of the specified order. If
154  // order < 1, we instead use an isoparametric/isogeometric space.
156  bool delete_fec;
157  if (order > 0)
158  {
159  fec = new H1_FECollection(order, dim);
160  delete_fec = true;
161  }
162  else if (pmesh.GetNodes())
163  {
164  fec = pmesh.GetNodes()->OwnFEC();
165  delete_fec = false;
166  if (myid == 0)
167  {
168  cout << "Using isoparametric FEs: " << fec->Name() << endl;
169  }
170  }
171  else
172  {
173  fec = new H1_FECollection(order = 1, dim);
174  delete_fec = true;
175  }
176  ParFiniteElementSpace fespace(&pmesh, fec);
177  HYPRE_BigInt size = fespace.GlobalTrueVSize();
178  if (myid == 0)
179  {
180  cout << "Number of finite element unknowns: " << size << endl;
181  }
182 
183  // 8. Determine the list of true (i.e. parallel conforming) essential
184  // boundary dofs. In this example, the boundary conditions are defined
185  // by marking all the boundary attributes from the mesh as essential
186  // (Dirichlet) and converting them to a list of true dofs.
187  Array<int> ess_tdof_list;
188  if (pmesh.bdr_attributes.Size())
189  {
190  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
191  ess_bdr = 1;
192  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
193  }
194 
195  // 9. Set up the parallel linear form b(.) which corresponds to the
196  // right-hand side of the FEM linear system, which in this case is
197  // (1,phi_i) where phi_i are the basis functions in fespace.
198  MFEM_PERF_BEGIN("Set up the linear form");
199  ParLinearForm b(&fespace);
200  ConstantCoefficient one(1.0);
202  b.Assemble();
203  MFEM_PERF_END("Set up the linear form");
204 
205  // 10. Define the solution vector x as a parallel finite element grid function
206  // corresponding to fespace. Initialize x with initial guess of zero,
207  // which satisfies the boundary conditions.
208  ParGridFunction x(&fespace);
209  x = 0.0;
210 
211  // 11. Set up the parallel bilinear form a(.,.) on the finite element space
212  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
213  // domain integrator.
214  MFEM_PERF_BEGIN("Set up the bilinear form");
215  ParBilinearForm a(&fespace);
216  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
218 
219  // 12. Assemble the parallel bilinear form and the corresponding linear
220  // system, applying any necessary transformations such as: parallel
221  // assembly, eliminating boundary conditions, applying conforming
222  // constraints for non-conforming AMR, static condensation, etc.
223  if (static_cond) { a.EnableStaticCondensation(); }
224  a.Assemble();
225 
226  OperatorPtr A;
227  Vector B, X;
228  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
229  MFEM_PERF_END("Set up the bilinear form");
230  // 13. Solve the linear system A X = B.
231  // * With full assembly, use the BoomerAMG preconditioner from hypre.
232  // * With partial assembly, use Jacobi smoothing, for now.
233  {
234  MFEM_PERF_SCOPE("Solve A X=B");
235  Solver *prec = NULL;
236  if (pa)
237  {
238  if (UsesTensorBasis(fespace))
239  {
240  prec = new OperatorJacobiSmoother(a, ess_tdof_list);
241  }
242  }
243  else
244  {
245  prec = new HypreBoomerAMG;
246  }
247  CGSolver cg(MPI_COMM_WORLD);
248  cg.SetRelTol(1e-12);
249  cg.SetMaxIter(2000);
250  cg.SetPrintLevel(1);
251  if (prec) { cg.SetPreconditioner(*prec); }
252  cg.SetOperator(*A);
253  cg.Mult(B, X);
254  delete prec;
255  }
256  // 14. Recover the parallel grid function corresponding to X. This is the
257  // local finite element solution on each processor.
258  a.RecoverFEMSolution(X, b, x);
259 
260  // 15. Save the refined mesh and the solution in parallel. This output can
261  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
262  MFEM_PERF_BEGIN("Save the results");
263  {
264  ostringstream mesh_name, sol_name;
265  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
266  sol_name << "sol." << setfill('0') << setw(6) << myid;
267 
268  ofstream mesh_ofs(mesh_name.str().c_str());
269  mesh_ofs.precision(8);
270  pmesh.Print(mesh_ofs);
271 
272  ofstream sol_ofs(sol_name.str().c_str());
273  sol_ofs.precision(8);
274  x.Save(sol_ofs);
275  }
276  MFEM_PERF_END("Save the results");
277  // 16. Send the solution by socket to a GLVis server.
278  if (visualization)
279  {
280  char vishost[] = "localhost";
281  int visport = 19916;
282  socketstream sol_sock(vishost, visport);
283  sol_sock << "parallel " << num_procs << " " << myid << "\n";
284  sol_sock.precision(8);
285  sol_sock << "solution\n" << pmesh << x << flush;
286  }
287 
288  // 17. Free the used memory.
289  if (delete_fec)
290  {
291  delete fec;
292  }
293  // Flush output before MPI_finalize
294  mgr.flush();
295 
296  return 0;
297 }
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:97
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:78
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
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
bool UsesTensorBasis(const FiniteElementSpace &fes)
Definition: fespace.hpp:983
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1443
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
constexpr char vishost[]
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
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.
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
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
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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
HYPRE_Int HYPRE_BigInt
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
FDualNumber< tbase > log(const FDualNumber< tbase > &f)
log([dual number])
Definition: fdual.hpp:515
double a
Definition: lissajous.cpp:41
virtual const char * Name() const
Definition: fe_coll.hpp:61
int dim
Definition: ex24.cpp:53
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.
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:904
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7841
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:216
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4684
Base class for solvers.
Definition: operator.hpp:651
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150