MFEM  v4.3.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 // 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
17 // space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
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);
59  args.AddOption(&mesh_file, "-m", "--mesh",
60  "Mesh file to use.");
61  args.AddOption(&order, "-o", "--order",
62  "Finite element order (polynomial degree) or -1 for"
63  " isoparametric space.");
64  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
65  "--no-static-condensation", "Enable static condensation.");
66  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
67  "--no-partial-assembly", "Enable Partial Assembly.");
68  args.AddOption(&amgx_lib, "-amgx", "--amgx-lib", "-no-amgx",
69  "--no-amgx-lib", "Use AmgX in example.");
70  args.AddOption(&amgx_json_file, "--amgx-file", "--amgx-file",
71  "AMGX solver config file (overrides --amgx-solver, --amgx-verbose)");
72  args.AddOption(&amgx_mpi_teams, "--amgx-mpi-teams", "--amgx-mpi-teams",
73  "--amgx-mpi-gpu-exclusive", "--amgx-mpi-gpu-exclusive",
74  "Create MPI teams when using AmgX to load balance between ranks and GPUs.");
75  args.AddOption(&device_config, "-d", "--device",
76  "Device configuration string, see Device::Configure().");
77  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
78  "--no-visualization",
79  "Enable or disable GLVis visualization.");
80  args.AddOption(&ndevices, "-nd","--gpus-per-node-in-teams-mode",
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_BigInt 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;
253  amgx.ReadParameters(amgx_json_file, AmgXSolver::EXTERNAL);
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.SetConvergenceCheck(true);
268  amgx.Mult(B, X);
269 
270  // Release MPI communicators and resources created by AmgX
271  amgx.Finalize();
272  }
273  else
274  {
275  prec = new HypreBoomerAMG;
276 
277  CGSolver cg(MPI_COMM_WORLD);
278  cg.SetRelTol(1e-12);
279  cg.SetMaxIter(2000);
280  cg.SetPrintLevel(1);
281  if (prec) { cg.SetPreconditioner(*prec); }
282  cg.SetOperator(*A);
283  cg.Mult(B, X);
284  delete prec;
285  }
286 
287  // 14. Recover the parallel grid function corresponding to X. This is the
288  // local finite element solution on each processor.
289  a.RecoverFEMSolution(X, b, x);
290 
291  // 15. Save the refined mesh and the solution in parallel. This output can
292  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
293  {
294  ostringstream mesh_name, sol_name;
295  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
296  sol_name << "sol." << setfill('0') << setw(6) << myid;
297 
298  ofstream mesh_ofs(mesh_name.str().c_str());
299  mesh_ofs.precision(8);
300  pmesh.Print(mesh_ofs);
301 
302  ofstream sol_ofs(sol_name.str().c_str());
303  sol_ofs.precision(8);
304  x.Save(sol_ofs);
305  }
306 
307  // 16. Send the solution by socket to a GLVis server.
308  if (visualization)
309  {
310  char vishost[] = "localhost";
311  int visport = 19916;
312  socketstream sol_sock(vishost, visport);
313  sol_sock << "parallel " << num_procs << " " << myid << "\n";
314  sol_sock.precision(8);
315  sol_sock << "solution\n" << pmesh << x << flush;
316  }
317 
318  // 17. Free the used memory.
319  if (delete_fec)
320  {
321  delete fec;
322  }
323  MPI_Finalize();
324 
325  return 0;
326 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
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:790
Conjugate gradient method.
Definition: solvers.hpp:316
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:99
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:616
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:275
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:841
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
void InitMPITeams(const MPI_Comm &comm, const int nDevs)
Definition: amgxsolver.cpp:149
bool UsesTensorBasis(const FiniteElementSpace &fes)
Definition: fespace.hpp:957
void ReadParameters(const std::string config, CONFIG_SRC source)
Definition: amgxsolver.cpp:186
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: amgxsolver.cpp:897
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
void SetConvergenceCheck(bool setConvergenceCheck_=true)
Add a check for convergence after applying Mult.
Definition: amgxsolver.cpp:193
Class for parallel linear form.
Definition: plinearform.hpp:26
void SetPrintLevel(int print_lvl)
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:150
constexpr char vishost[]
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:130
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
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.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4382
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
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:204
void SetRelTol(double rtol)
Definition: solvers.hpp:98
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:854
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)
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:120
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:327
Class for parallel bilinear form.
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:828
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:330
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:92
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Base class for solvers.
Definition: operator.hpp:648
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
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
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