MFEM  v4.2.0 Finite element discretization library
ex1p.cpp
Go to the documentation of this file.
1 // MFEM Example 1 - Parallel Version
2 // AmgX Modification
3 //
4 // Compile with: make ex1p
5 //
6 // AmgX sample runs:
7 // mpirun -np 4 ex1p
8 // mpirun -np 4 ex1p -d cuda
9 // mpirun -np 10 ex1p --amgx-file amg_pcg.json --amgx-mpi-teams
10 // mpirun -np 4 ex1p --amgx-file amg_pcg.json
11 //
12 // Description: This example code demonstrates the use of MFEM to define a
13 // simple finite element discretization of the Laplace problem
14 // -Delta u = 1 with homogeneous Dirichlet boundary conditions.
15 // Specifically, we discretize using a FE space of the specified
16 // order, or if order < 1 using an isoparametric/isogeometric
18 // NURBS mesh, etc.)
19 //
20 // The example highlights the use of mesh refinement, finite
21 // element grid functions, as well as linear and bilinear forms
22 // corresponding to the left-hand side and right-hand side of the
23 // discrete linear system. We also cover the explicit elimination
24 // of essential boundary conditions, static condensation, and the
25 // optional connection to the GLVis tool for visualization.
26
27 #include "mfem.hpp"
28 #include <fstream>
29 #include <iostream>
30
31 using namespace std;
32 using namespace mfem;
33
34 #ifndef MFEM_USE_AMGX
35 #error This example requires that MFEM is built with MFEM_USE_AMGX=YES
36 #endif
37
38 int main(int argc, char *argv[])
39 {
40  // 1. Initialize MPI.
41  int num_procs, myid;
42  MPI_Init(&argc, &argv);
43  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
44  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
45
46  // 2. Parse command-line options.
47  const char *mesh_file = "../../data/star.mesh";
48  int order = 1;
49  bool static_cond = false;
50  bool pa = false;
51  const char *device_config = "cpu";
52  bool visualization = true;
53  bool amgx_lib = true;
54  bool amgx_mpi_teams = false;
55  const char* amgx_json_file = ""; // JSON file for AmgX
56  int ndevices = 1;
57
58  OptionsParser args(argc, argv);
60  "Mesh file to use.");
62  "Finite element order (polynomial degree) or -1 for"
63  " isoparametric space.");
65  "--no-static-condensation", "Enable static condensation.");
67  "--no-partial-assembly", "Enable Partial Assembly.");
69  "--no-amgx-lib", "Use AmgX in example.");
71  "AMGX solver config file (overrides --amgx-solver, --amgx-verbose)");
73  "--amgx-mpi-gpu-exclusive", "--amgx-mpi-gpu-exclusive",
74  "Create MPI teams when using AmgX to load balance between ranks and GPUs.");
76  "Device configuration string, see Device::Configure().");
78  "--no-visualization",
79  "Enable or disable GLVis visualization.");
81  "Number of GPU devices per node (Only used if amgx_mpi_teams is true).");
82
83  args.Parse();
84  if (!args.Good())
85  {
86  if (myid == 0)
87  {
88  args.PrintUsage(cout);
89  }
90  MPI_Finalize();
91  return 1;
92  }
93  if (myid == 0)
94  {
95  args.PrintOptions(cout);
96  }
97
98  // 3. Enable hardware devices such as GPUs, and programming models such as
99  // CUDA, OCCA, RAJA and OpenMP based on command line options.
100  Device device(device_config);
101  if (myid == 0) { device.Print(); }
102
103  // 4. Read the (serial) mesh from the given mesh file on all processors. We
104  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
105  // and volume meshes with the same code.
106  Mesh mesh(mesh_file, 1, 1);
107  int dim = mesh.Dimension();
108
109  // 5. Refine the serial mesh on all processors to increase the resolution. In
110  // this example we do 'ref_levels' of uniform refinement. We choose
111  // 'ref_levels' to be the largest number that gives a final mesh with no
112  // more than 10,000 elements.
113  {
114  int ref_levels =
115  (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
116  for (int l = 0; l < ref_levels; l++)
117  {
118  mesh.UniformRefinement();
119  }
120  }
121
122  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
123  // this mesh further in parallel to increase the resolution. Once the
124  // parallel mesh is defined, the serial mesh can be deleted.
125  ParMesh pmesh(MPI_COMM_WORLD, mesh);
126  mesh.Clear();
127  {
128  int par_ref_levels = 2;
129  for (int l = 0; l < par_ref_levels; l++)
130  {
131  pmesh.UniformRefinement();
132  }
133  }
134
135  // 7. Define a parallel finite element space on the parallel mesh. Here we
136  // use continuous Lagrange finite elements of the specified order. If
137  // order < 1, we instead use an isoparametric/isogeometric space.
139  bool delete_fec;
140  if (order > 0)
141  {
142  fec = new H1_FECollection(order, dim);
143  delete_fec = true;
144  }
145  else if (pmesh.GetNodes())
146  {
147  fec = pmesh.GetNodes()->OwnFEC();
148  delete_fec = false;
149  if (myid == 0)
150  {
151  cout << "Using isoparametric FEs: " << fec->Name() << endl;
152  }
153  }
154  else
155  {
156  fec = new H1_FECollection(order = 1, dim);
157  delete_fec = true;
158  }
159  ParFiniteElementSpace fespace(&pmesh, fec);
160  HYPRE_Int size = fespace.GlobalTrueVSize();
161  if (myid == 0)
162  {
163  cout << "Number of finite element unknowns: " << size << endl;
164  }
165
166  // 8. Determine the list of true (i.e. parallel conforming) essential
167  // boundary dofs. In this example, the boundary conditions are defined
168  // by marking all the boundary attributes from the mesh as essential
169  // (Dirichlet) and converting them to a list of true dofs.
170  Array<int> ess_tdof_list;
171  if (pmesh.bdr_attributes.Size())
172  {
173  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
174  ess_bdr = 1;
175  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
176  }
177
178  // 9. Set up the parallel linear form b(.) which corresponds to the
179  // right-hand side of the FEM linear system, which in this case is
180  // (1,phi_i) where phi_i are the basis functions in fespace.
181  ParLinearForm b(&fespace);
182  ConstantCoefficient one(1.0);
184  b.Assemble();
185
186  // 10. Define the solution vector x as a parallel finite element grid function
187  // corresponding to fespace. Initialize x with initial guess of zero,
188  // which satisfies the boundary conditions.
189  ParGridFunction x(&fespace);
190  x = 0.0;
191
192  // 11. Set up the parallel bilinear form a(.,.) on the finite element space
193  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
194  // domain integrator.
195  ParBilinearForm a(&fespace);
196  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
198
199  // 12. Assemble the parallel bilinear form and the corresponding linear
200  // system, applying any necessary transformations such as: parallel
201  // assembly, eliminating boundary conditions, applying conforming
202  // constraints for non-conforming AMR, static condensation, etc.
203  if (static_cond) { a.EnableStaticCondensation(); }
204  a.Assemble();
205
206  OperatorPtr A;
207  Vector B, X;
208  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
209
210  // 13. Solve the linear system A X = B.
211  // * With full assembly, use the BoomerAMG preconditioner from hypre.
212  // * If AmgX is available solve using amg preconditioner.
213  // * With partial assembly, use Jacobi smoothing, for now.
214  Solver *prec = NULL;
215  if (pa)
216  {
217  if (UsesTensorBasis(fespace))
218  {
219  prec = new OperatorJacobiSmoother(a, ess_tdof_list);
220  }
221
222  CGSolver cg(MPI_COMM_WORLD);
223  cg.SetRelTol(1e-12);
224  cg.SetMaxIter(2000);
225  cg.SetPrintLevel(1);
226  if (prec) { cg.SetPreconditioner(*prec); }
227  cg.SetOperator(*A);
228  cg.Mult(B, X);
229  delete prec;
230  }
231  else if (amgx_lib && strcmp(amgx_json_file,"") == 0)
232  {
233  MFEM_VERIFY(!amgx_mpi_teams,
234  "Please add JSON file to try AmgX with MPI teams mode");
235
236  bool amgx_verbose = false;
237  prec = new AmgXSolver(MPI_COMM_WORLD, AmgXSolver::PRECONDITIONER,
238  amgx_verbose);
239
240  CGSolver cg(MPI_COMM_WORLD);
241  cg.SetRelTol(1e-12);
242  cg.SetMaxIter(2000);
243  cg.SetPrintLevel(1);
244  if (prec) { cg.SetPreconditioner(*prec); }
245  cg.SetOperator(*A);
246  cg.Mult(B, X);
247  delete prec;
248
249  }
250  else if (amgx_lib && strcmp(amgx_json_file,"") != 0)
251  {
252  AmgXSolver amgx;
254
255  if (amgx_mpi_teams)
256  {
257  // Forms MPI teams to load balance between MPI ranks and GPUs
258  amgx.InitMPITeams(MPI_COMM_WORLD, ndevices);
259  }
260  else
261  {
262  // Assumes each MPI rank is paired with a GPU
263  amgx.InitExclusiveGPU(MPI_COMM_WORLD);
264  }
265
266  amgx.SetOperator(*A.As<HypreParMatrix>());
267  amgx.Mult(B, X);
268
269  // Release MPI communicators and resources created by AmgX
270  amgx.Finalize();
271  }
272  else
273  {
274  prec = new HypreBoomerAMG;
275
276  CGSolver cg(MPI_COMM_WORLD);
277  cg.SetRelTol(1e-12);
278  cg.SetMaxIter(2000);
279  cg.SetPrintLevel(1);
280  if (prec) { cg.SetPreconditioner(*prec); }
281  cg.SetOperator(*A);
282  cg.Mult(B, X);
283  delete prec;
284  }
285
286  // 14. Recover the parallel grid function corresponding to X. This is the
287  // local finite element solution on each processor.
288  a.RecoverFEMSolution(X, b, x);
289
290  // 15. Save the refined mesh and the solution in parallel. This output can
291  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
292  {
293  ostringstream mesh_name, sol_name;
294  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
295  sol_name << "sol." << setfill('0') << setw(6) << myid;
296
297  ofstream mesh_ofs(mesh_name.str().c_str());
298  mesh_ofs.precision(8);
299  pmesh.Print(mesh_ofs);
300
301  ofstream sol_ofs(sol_name.str().c_str());
302  sol_ofs.precision(8);
303  x.Save(sol_ofs);
304  }
305
306  // 16. Send the solution by socket to a GLVis server.
307  if (visualization)
308  {
309  char vishost[] = "localhost";
310  int visport = 19916;
311  socketstream sol_sock(vishost, visport);
312  sol_sock << "parallel " << num_procs << " " << myid << "\n";
313  sol_sock.precision(8);
314  sol_sock << "solution\n" << pmesh << x << flush;
315  }
316
317  // 17. Free the used memory.
318  if (delete_fec)
319  {
320  delete fec;
321  }
322  MPI_Finalize();
323
324  return 0;
325 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:93
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:775
Definition: solvers.hpp:258
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:96
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:535
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:737
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:819
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:261
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void InitMPITeams(const MPI_Comm &comm, const int nDevs)
Definition: amgxsolver.cpp:137
bool UsesTensorBasis(const FiniteElementSpace &fes)
Definition: fespace.hpp:983
int main(int argc, char *argv[])
Definition: ex1.cpp:66
void ReadParameters(const std::string config, CONFIG_SRC source)
Definition: amgxsolver.cpp:174
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: amgxsolver.cpp:864
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1079
Class for parallel linear form.
Definition: plinearform.hpp:26
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:70
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
constexpr char vishost[]
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:128
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:98
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:251
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4169
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
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:201
void SetRelTol(double rtol)
Definition: solvers.hpp:96
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
virtual void SetOperator(const Operator &op)
Definition: amgxsolver.cpp:821
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
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
double a
Definition: lissajous.cpp:41
virtual const char * Name() const
Definition: fe_coll.hpp:61
void InitExclusiveGPU(const MPI_Comm &comm)
Definition: amgxsolver.cpp:108
int dim
Definition: ex24.cpp:53
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:304
Class for parallel bilinear form.
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:719
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:272
Vector data type.
Definition: vector.hpp:51
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6603
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:91
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:159
Base class for solvers.
Definition: operator.hpp:634
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:118
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:181
Class for parallel meshes.
Definition: pmesh.hpp:32
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:145