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 // 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 and HYPRE.
41  Mpi::Init(argc, argv);
42  int num_procs = Mpi::WorldSize();
43  int myid = Mpi::WorldRank();
44  Hypre::Init();
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  return 1;
91  }
92  if (myid == 0)
93  {
94  args.PrintOptions(cout);
95  }
96 
97  // 3. Enable hardware devices such as GPUs, and programming models such as
98  // CUDA, OCCA, RAJA and OpenMP based on command line options.
99  Device device(device_config);
100  if (myid == 0) { device.Print(); }
101 
102  // 4. Read the (serial) mesh from the given mesh file on all processors. We
103  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
104  // and volume meshes with the same code.
105  Mesh mesh(mesh_file, 1, 1);
106  int dim = mesh.Dimension();
107 
108  // 5. Refine the serial mesh on all processors to increase the resolution. In
109  // this example we do 'ref_levels' of uniform refinement. We choose
110  // 'ref_levels' to be the largest number that gives a final mesh with no
111  // more than 10,000 elements.
112  {
113  int ref_levels =
114  (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
115  for (int l = 0; l < ref_levels; l++)
116  {
117  mesh.UniformRefinement();
118  }
119  }
120 
121  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
122  // this mesh further in parallel to increase the resolution. Once the
123  // parallel mesh is defined, the serial mesh can be deleted.
124  ParMesh pmesh(MPI_COMM_WORLD, mesh);
125  mesh.Clear();
126  {
127  int par_ref_levels = 2;
128  for (int l = 0; l < par_ref_levels; l++)
129  {
130  pmesh.UniformRefinement();
131  }
132  }
133 
134  // 7. Define a parallel finite element space on the parallel mesh. Here we
135  // use continuous Lagrange finite elements of the specified order. If
136  // order < 1, we instead use an isoparametric/isogeometric space.
138  bool delete_fec;
139  if (order > 0)
140  {
141  fec = new H1_FECollection(order, dim);
142  delete_fec = true;
143  }
144  else if (pmesh.GetNodes())
145  {
146  fec = pmesh.GetNodes()->OwnFEC();
147  delete_fec = false;
148  if (myid == 0)
149  {
150  cout << "Using isoparametric FEs: " << fec->Name() << endl;
151  }
152  }
153  else
154  {
155  fec = new H1_FECollection(order = 1, dim);
156  delete_fec = true;
157  }
158  ParFiniteElementSpace fespace(&pmesh, fec);
159  HYPRE_BigInt size = fespace.GlobalTrueVSize();
160  if (myid == 0)
161  {
162  cout << "Number of finite element unknowns: " << size << endl;
163  }
164 
165  // 8. Determine the list of true (i.e. parallel conforming) essential
166  // boundary dofs. In this example, the boundary conditions are defined
167  // by marking all the boundary attributes from the mesh as essential
168  // (Dirichlet) and converting them to a list of true dofs.
169  Array<int> ess_tdof_list;
170  if (pmesh.bdr_attributes.Size())
171  {
172  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
173  ess_bdr = 1;
174  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
175  }
176 
177  // 9. Set up the parallel linear form b(.) which corresponds to the
178  // right-hand side of the FEM linear system, which in this case is
179  // (1,phi_i) where phi_i are the basis functions in fespace.
180  ParLinearForm b(&fespace);
181  ConstantCoefficient one(1.0);
183  b.Assemble();
184 
185  // 10. Define the solution vector x as a parallel finite element grid function
186  // corresponding to fespace. Initialize x with initial guess of zero,
187  // which satisfies the boundary conditions.
188  ParGridFunction x(&fespace);
189  x = 0.0;
190 
191  // 11. Set up the parallel bilinear form a(.,.) on the finite element space
192  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
193  // domain integrator.
194  ParBilinearForm a(&fespace);
195  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
197 
198  // 12. Assemble the parallel bilinear form and the corresponding linear
199  // system, applying any necessary transformations such as: parallel
200  // assembly, eliminating boundary conditions, applying conforming
201  // constraints for non-conforming AMR, static condensation, etc.
202  if (static_cond) { a.EnableStaticCondensation(); }
203  a.Assemble();
204 
205  OperatorPtr A;
206  Vector B, X;
207  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
208 
209  // 13. Solve the linear system A X = B.
210  // * With full assembly, use the BoomerAMG preconditioner from hypre.
211  // * If AmgX is available solve using amg preconditioner.
212  // * With partial assembly, use Jacobi smoothing, for now.
213  Solver *prec = NULL;
214  if (pa)
215  {
216  if (UsesTensorBasis(fespace))
217  {
218  prec = new OperatorJacobiSmoother(a, ess_tdof_list);
219  }
220 
221  CGSolver cg(MPI_COMM_WORLD);
222  cg.SetRelTol(1e-12);
223  cg.SetMaxIter(2000);
224  cg.SetPrintLevel(1);
225  if (prec) { cg.SetPreconditioner(*prec); }
226  cg.SetOperator(*A);
227  cg.Mult(B, X);
228  delete prec;
229  }
230  else if (amgx_lib && strcmp(amgx_json_file,"") == 0)
231  {
232  MFEM_VERIFY(!amgx_mpi_teams,
233  "Please add JSON file to try AmgX with MPI teams mode");
234 
235  bool amgx_verbose = false;
236  prec = new AmgXSolver(MPI_COMM_WORLD, AmgXSolver::PRECONDITIONER,
237  amgx_verbose);
238 
239  CGSolver cg(MPI_COMM_WORLD);
240  cg.SetRelTol(1e-12);
241  cg.SetMaxIter(2000);
242  cg.SetPrintLevel(1);
243  if (prec) { cg.SetPreconditioner(*prec); }
244  cg.SetOperator(*A);
245  cg.Mult(B, X);
246  delete prec;
247 
248  }
249  else if (amgx_lib && strcmp(amgx_json_file,"") != 0)
250  {
251  AmgXSolver amgx;
252  amgx.ReadParameters(amgx_json_file, AmgXSolver::EXTERNAL);
253 
254  if (amgx_mpi_teams)
255  {
256  // Forms MPI teams to load balance between MPI ranks and GPUs
257  amgx.InitMPITeams(MPI_COMM_WORLD, ndevices);
258  }
259  else
260  {
261  // Assumes each MPI rank is paired with a GPU
262  amgx.InitExclusiveGPU(MPI_COMM_WORLD);
263  }
264 
265  amgx.SetOperator(*A.As<HypreParMatrix>());
266  amgx.SetConvergenceCheck(true);
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 
323  return 0;
324 }
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
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
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
void InitMPITeams(const MPI_Comm &comm, const int nDevs)
Definition: amgxsolver.cpp:149
bool UsesTensorBasis(const FiniteElementSpace &fes)
Definition: fespace.hpp:983
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:902
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1443
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
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
virtual void SetOperator(const Operator &op)
Definition: amgxsolver.cpp:859
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
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: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
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:337
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