MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex1.cpp
Go to the documentation of this file.
1// MFEM Example 1
2// AmgX Modification
3//
4// Compile with: make ex1
5//
6// AmgX sample runs:
7// ex1
8// ex1 -d cuda
9// ex1 --amgx-file multi_gs.json --amgx-solver
10// ex1 --amgx-file precon.json --amgx-preconditioner
11// ex1 --amgx-file multi_gs.json --amgx-solver -d cuda
12// ex1 --amgx-file precon.json --amgx-preconditioner -d cuda
13//
14// Description: This example code demonstrates the use of MFEM to define a
15// simple finite element discretization of the Laplace problem
16// -Delta u = 1 with homogeneous Dirichlet boundary conditions.
17// Specifically, we discretize using a FE space of the specified
18// order, or if order < 1 using an isoparametric/isogeometric
20// NURBS mesh, etc.)
21//
22// The example highlights the use of mesh refinement, finite
23// element grid functions, as well as linear and bilinear forms
24// corresponding to the left-hand side and right-hand side of the
25// discrete linear system. We also cover the explicit elimination
26// of essential boundary conditions, static condensation, and the
27// optional connection to the GLVis tool for visualization.
28
29#include "mfem.hpp"
30#include <fstream>
31#include <iostream>
32
33#ifndef MFEM_USE_AMGX
34#error This example requires that MFEM is built with MFEM_USE_AMGX=YES
35#endif
36
37using namespace std;
38using namespace mfem;
39
40int main(int argc, char *argv[])
41{
42 // 1. Parse command-line options.
43 const char *mesh_file = "../../data/star.mesh";
44 int order = 1;
45 bool static_cond = false;
46 bool pa = false;
47 const char *device_config = "cpu";
48 bool visualization = true;
49 bool amgx_lib = true;
50 bool amgx_solver = true;
51 const char* amgx_json_file = ""; // JSON file for AmgX
52
53 OptionsParser args(argc, argv);
55 "Mesh file to use.");
57 "Finite element order (polynomial degree) or -1 for"
58 " isoparametric space.");
60 "--no-static-condensation", "Enable static condensation.");
62 "--no-partial-assembly", "Enable Partial Assembly.");
64 "--no-amgx-lib", "Use AmgX in example.");
66 "AMGX solver config file (overrides --amgx-solver, --amgx-verbose)");
68 "--amgx-preconditioner", "--amgx-preconditioner",
69 "Configure AMGX as solver or preconditioner.");
71 "Device configuration string, see Device::Configure().");
73 "--no-visualization",
74 "Enable or disable GLVis visualization.");
75 args.Parse();
76 if (!args.Good())
77 {
78 args.PrintUsage(cout);
79 return 1;
80 }
81 args.PrintOptions(cout);
82
83 // 2. Enable hardware devices such as GPUs, and programming models such as
84 // CUDA, OCCA, RAJA and OpenMP based on command line options.
85 Device device(device_config);
86 device.Print();
87
88 // 3. Read the mesh from the given mesh file. We can handle triangular,
89 // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
90 // the same code.
91 Mesh mesh(mesh_file, 1, 1);
92 int dim = mesh.Dimension();
93
94 // 4. Refine the mesh to increase the resolution. In this example we do
95 // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
96 // largest number that gives a final mesh with no more than 50,000
97 // elements.
98 {
99 int ref_levels =
100 (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
101 for (int l = 0; l < ref_levels; l++)
102 {
103 mesh.UniformRefinement();
104 }
105 }
106
107 // 5. Define a finite element space on the mesh. Here we use continuous
108 // Lagrange finite elements of the specified order. If order < 1, we
109 // instead use an isoparametric/isogeometric space.
111 bool delete_fec;
112 if (order > 0)
113 {
114 fec = new H1_FECollection(order, dim);
115 delete_fec = true;
116 }
117 else if (mesh.GetNodes())
118 {
119 fec = mesh.GetNodes()->OwnFEC();
120 delete_fec = false;
121 cout << "Using isoparametric FEs: " << fec->Name() << endl;
122 }
123 else
124 {
125 fec = new H1_FECollection(order = 1, dim);
126 delete_fec = true;
127 }
128 FiniteElementSpace fespace(&mesh, fec);
129 cout << "Number of finite element unknowns: "
130 << fespace.GetTrueVSize() << endl;
131
132 // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
133 // In this example, the boundary conditions are defined by marking all
134 // the boundary attributes from the mesh as essential (Dirichlet) and
135 // converting them to a list of true dofs.
136 Array<int> ess_tdof_list;
137 if (mesh.bdr_attributes.Size())
138 {
139 Array<int> ess_bdr(mesh.bdr_attributes.Max());
140 ess_bdr = 1;
141 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
142 }
143
144 // 7. Set up the linear form b(.) which corresponds to the right-hand side of
145 // the FEM linear system, which in this case is (1,phi_i) where phi_i are
146 // the basis functions in the finite element fespace.
147 LinearForm b(&fespace);
148 ConstantCoefficient one(1.0);
150 b.Assemble();
151
152 // 8. Define the solution vector x as a finite element grid function
153 // corresponding to fespace. Initialize x with initial guess of zero,
154 // which satisfies the boundary conditions.
155 GridFunction x(&fespace);
156 x = 0.0;
157
158 // 9. Set up the bilinear form a(.,.) on the finite element space
159 // corresponding to the Laplacian operator -Delta, by adding the Diffusion
160 // domain integrator.
161 BilinearForm a(&fespace);
162 if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
164
165 // 10. Assemble the bilinear form and the corresponding linear system,
166 // applying any necessary transformations such as: eliminating boundary
167 // conditions, applying conforming constraints for non-conforming AMR,
168 // static condensation, etc.
169 if (static_cond) { a.EnableStaticCondensation(); }
170 a.Assemble();
171
172 OperatorPtr A;
173 Vector B, X;
174 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
175
176 cout << "Size of linear system: " << A->Height() << endl;
177
178 // 11. Solve the linear system A X = B.
179 if (pa)
180 {
181 // Jacobi preconditioning in partial assembly mode
182 if (UsesTensorBasis(fespace))
183 {
184 OperatorJacobiSmoother M(a, ess_tdof_list);
185 PCG(*A, M, B, X, 1, 400, 1e-12, 0.0);
186 }
187 else
188 {
189 CG(*A, B, X, 1, 400, 1e-12, 0.0);
190 }
191 }
192 else if (amgx_lib && strcmp(amgx_json_file,"") == 0)
193 {
194 bool amgx_verbose = false;
195 AmgXSolver amgx(AmgXSolver::PRECONDITIONER, amgx_verbose);
196 amgx.SetOperator(*A.As<SparseMatrix>());
197 PCG(*A, amgx, B, X, 1, 200, 1e-12, 0.0);
198 }
199 else if (amgx_lib && strcmp(amgx_json_file,"") != 0)
200 {
201 AmgXSolver amgx;
203 amgx.InitSerial();
204 amgx.SetOperator(*A.As<SparseMatrix>());
205
206 if (amgx_solver)
207 {
208 amgx.SetConvergenceCheck(true);
209 amgx.Mult(B,X);
210 }
211 else
212 {
213 // Omit convergence check at the AmgX level when using as a
214 // preconditioner.
215 amgx.SetConvergenceCheck(false);
216 PCG(*A.As<SparseMatrix>(), amgx, B, X, 3, 40, 1e-12, 0.0);
217 }
218 }
219 else
220 {
221#ifndef MFEM_USE_SUITESPARSE
222 // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
223 GSSmoother M((SparseMatrix&)(*A));
224 PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
225#else
226 // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
227 UMFPackSolver umf_solver;
228 umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
229 umf_solver.SetOperator(*A);
230 umf_solver.Mult(B, X);
231#endif
232 }
233
234 // 12. Recover the solution as a finite element grid function.
235 a.RecoverFEMSolution(X, b, x);
236
237 // 13. Save the refined mesh and the solution. This output can be viewed later
238 // using GLVis: "glvis -m refined.mesh -g sol.gf".
239 ofstream mesh_ofs("refined.mesh");
240 mesh_ofs.precision(8);
241 mesh.Print(mesh_ofs);
242 ofstream sol_ofs("sol.gf");
243 sol_ofs.precision(8);
244 x.Save(sol_ofs);
245
246 // 14. Send the solution by socket to a GLVis server.
247 if (visualization)
248 {
249 char vishost[] = "localhost";
250 int visport = 19916;
251 socketstream sol_sock(vishost, visport);
252 sol_sock.precision(8);
253 sol_sock << "solution\n" << mesh << x << flush;
254 }
255
256 // 15. Free the used memory.
257 if (delete_fec)
258 {
259 delete fec;
260 }
261
262 return 0;
263}
MFEM wrapper for Nvidia's multigrid library, AmgX (github.com/NVIDIA/AMGX)
@ EXTERNAL
Configure will be read from a specified file.
virtual void SetOperator(const Operator &op)
Sets the Operator that is going to be solved via AmgX. Supports operators based on either an MFEM Spa...
void InitSerial()
Initialize the AmgX library for serial execution once the solver configuration has been established t...
void SetConvergenceCheck(bool setConvergenceCheck_=true)
Add a check for convergence after applying Mult.
virtual void Mult(const Vector &b, Vector &x) const
Utilize the AmgX library to solve the linear system where the "matrix" is the AMG approximation to th...
void ReadParameters(const std::string config, CONFIG_SRC source)
Read in the AmgX parameters either through a file or directly through a properly formated string....
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Size() const
Return the logical size of the array.
Definition array.hpp:144
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
A coefficient that is constant across space and time.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
Class for domain integration .
Definition lininteg.hpp:109
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
virtual const char * Name() const
Definition fe_coll.hpp:79
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:588
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Vector with associated FE space and LinearFormIntegrators.
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:313
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Data type sparse matrix.
Definition sparsemat.hpp:51
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1096
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition solvers.cpp:3102
real_t Control[UMFPACK_CONTROL]
Definition solvers.hpp:1106
virtual void Mult(const Vector &b, Vector &x) const
Direct solution of the linear system using UMFPACK.
Definition solvers.cpp:3197
Vector data type.
Definition vector.hpp:80
int dim
Definition ex24.cpp:53
int main()
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:913
void CG(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:898
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition fespace.hpp:1345
const char vishost[]