MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
ex8p.cpp
Go to the documentation of this file.
1// MFEM Example 8 - Parallel Version
2//
3// Compile with: make ex8p
4//
5// Sample runs: mpirun -np 4 ex8p -m ../data/square-disc.mesh
6// mpirun -np 4 ex8p -m ../data/star.mesh
7// mpirun -np 4 ex8p -m ../data/star-mixed.mesh
8// mpirun -np 4 ex8p -m ../data/escher.mesh
9// mpirun -np 4 ex8p -m ../data/fichera.mesh
10// mpirun -np 4 ex8p -m ../data/fichera-mixed.mesh
11// mpirun -np 4 ex8p -m ../data/square-disc-p2.vtk
12// mpirun -np 4 ex8p -m ../data/square-disc-p3.mesh
13// mpirun -np 4 ex8p -m ../data/star-surf.mesh -o 2
14//
15// Description: This example code demonstrates the use of the Discontinuous
16// Petrov-Galerkin (DPG) method in its primal 2x2 block form as a
17// simple finite element discretization of the Laplace problem
18// -Delta u = f with homogeneous Dirichlet boundary conditions. We
19// use high-order continuous trial space, a high-order interfacial
20// (trace) space, and a high-order discontinuous test space
21// defining a local dual (H^{-1}) norm.
22//
23// We use the primal form of DPG, see "A primal DPG method without
24// a first-order reformulation", Demkowicz and Gopalakrishnan, CAM
25// 2013, DOI:10.1016/j.camwa.2013.06.029.
26//
27// The example highlights the use of interfacial (trace) finite
28// elements and spaces, trace face integrators and the definition
29// of block operators and preconditioners. The use of the ADS
30// preconditioner from hypre for interfacially-reduced H(div)
31// problems is also illustrated.
32//
33// We recommend viewing examples 1-5 before viewing this example.
34
35#include "mfem.hpp"
36#include <fstream>
37#include <iostream>
38
39using namespace std;
40using namespace mfem;
41
42int main(int argc, char *argv[])
43{
44 // 1. Initialize MPI and HYPRE.
45 Mpi::Init(argc, argv);
46 int num_procs = Mpi::WorldSize();
47 int myid = Mpi::WorldRank();
49
50 // 2. Parse command-line options.
51 const char *mesh_file = "../data/star.mesh";
52 int order = 1;
53 bool visualization = 1;
54
55 OptionsParser args(argc, argv);
56 args.AddOption(&mesh_file, "-m", "--mesh",
57 "Mesh file to use.");
58 args.AddOption(&order, "-o", "--order",
59 "Finite element order (polynomial degree).");
60 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
61 "--no-visualization",
62 "Enable or disable GLVis visualization.");
63 args.Parse();
64 if (!args.Good())
65 {
66 if (myid == 0)
67 {
68 args.PrintUsage(cout);
69 }
70 return 1;
71 }
72 if (myid == 0)
73 {
74 args.PrintOptions(cout);
75 }
76
77 // 3. Read the (serial) mesh from the given mesh file on all processors. We
78 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
79 // and volume meshes with the same code.
80 Mesh *mesh = new Mesh(mesh_file, 1, 1);
81 int dim = mesh->Dimension();
82
83 // 4. Refine the serial mesh on all processors to increase the resolution. In
84 // this example we do 'ref_levels' of uniform refinement. We choose
85 // 'ref_levels' to be the largest number that gives a final mesh with no
86 // more than 10,000 elements.
87 {
88 int ref_levels =
89 (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
90 for (int l = 0; l < ref_levels; l++)
91 {
92 mesh->UniformRefinement();
93 }
94 }
95
96 // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
97 // this mesh further in parallel to increase the resolution. Once the
98 // parallel mesh is defined, the serial mesh can be deleted.
99 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
100 delete mesh;
101 {
102 int par_ref_levels = 1;
103 for (int l = 0; l < par_ref_levels; l++)
104 {
105 pmesh->UniformRefinement();
106 }
107 }
108
109 // 6. Define the trial, interfacial (trace) and test DPG spaces:
110 // - The trial space, x0_space, contains the non-interfacial unknowns and
111 // has the essential BC.
112 // - The interfacial space, xhat_space, contains the interfacial unknowns
113 // and does not have essential BC.
114 // - The test space, test_space, is an enriched space where the enrichment
115 // degree may depend on the spatial dimension of the domain, the type of
116 // the mesh and the trial space order.
117 unsigned int trial_order = order;
118 unsigned int trace_order = order - 1;
119 unsigned int test_order = order; /* reduced order, full order is
120 (order + dim - 1) */
121 if (dim == 2 && (order%2 == 0 || (pmesh->MeshGenerator() & 2 && order > 1)))
122 {
123 test_order++;
124 }
125 if (test_order < trial_order)
126 {
127 if (myid == 0)
128 {
129 cerr << "Warning, test space not enriched enough to handle primal"
130 << " trial space\n";
131 }
132 }
133
134 FiniteElementCollection *x0_fec, *xhat_fec, *test_fec;
135
136 x0_fec = new H1_FECollection(trial_order, dim);
137 xhat_fec = new RT_Trace_FECollection(trace_order, dim);
138 test_fec = new L2_FECollection(test_order, dim);
139
140 ParFiniteElementSpace *x0_space, *xhat_space, *test_space;
141
142 x0_space = new ParFiniteElementSpace(pmesh, x0_fec);
143 xhat_space = new ParFiniteElementSpace(pmesh, xhat_fec);
144 test_space = new ParFiniteElementSpace(pmesh, test_fec);
145
146 HYPRE_BigInt glob_true_s0 = x0_space->GlobalTrueVSize();
147 HYPRE_BigInt glob_true_s1 = xhat_space->GlobalTrueVSize();
148 HYPRE_BigInt glob_true_s_test = test_space->GlobalTrueVSize();
149 if (myid == 0)
150 {
151 cout << "\nNumber of Unknowns:\n"
152 << " Trial space, X0 : " << glob_true_s0
153 << " (order " << trial_order << ")\n"
154 << " Interface space, Xhat : " << glob_true_s1
155 << " (order " << trace_order << ")\n"
156 << " Test space, Y : " << glob_true_s_test
157 << " (order " << test_order << ")\n\n";
158 }
159
160 // 7. Set up the linear form F(.) which corresponds to the right-hand side of
161 // the FEM linear system, which in this case is (f,phi_i) where f=1.0 and
162 // phi_i are the basis functions in the test finite element fespace.
163 ConstantCoefficient one(1.0);
164 ParLinearForm * F = new ParLinearForm(test_space);
166 F->Assemble();
167
168 ParGridFunction * x0 = new ParGridFunction(x0_space);
169 *x0 = 0.;
170
171 // 8. Set up the mixed bilinear form for the primal trial unknowns, B0,
172 // the mixed bilinear form for the interfacial unknowns, Bhat,
173 // the inverse stiffness matrix on the discontinuous test space, Sinv,
174 // and the stiffness matrix on the continuous trial space, S0.
175 Array<int> ess_bdr(pmesh->bdr_attributes.Max());
176 ess_bdr = 1;
177 Array<int> ess_dof;
178 x0_space->GetEssentialVDofs(ess_bdr, ess_dof);
179
180 ParMixedBilinearForm *B0 = new ParMixedBilinearForm(x0_space,test_space);
182 B0->Assemble();
183 B0->EliminateEssentialBCFromTrialDofs(ess_dof, *x0, *F);
184 B0->Finalize();
185
186 ParMixedBilinearForm *Bhat = new ParMixedBilinearForm(xhat_space,test_space);
188 Bhat->Assemble();
189 Bhat->Finalize();
190
191 ParBilinearForm *Sinv = new ParBilinearForm(test_space);
192 SumIntegrator *Sum = new SumIntegrator;
193 Sum->AddIntegrator(new DiffusionIntegrator(one));
194 Sum->AddIntegrator(new MassIntegrator(one));
196 Sinv->Assemble();
197 Sinv->Finalize();
198
199 ParBilinearForm *S0 = new ParBilinearForm(x0_space);
201 S0->Assemble();
202 S0->EliminateEssentialBC(ess_bdr);
203 S0->Finalize();
204
205 HypreParMatrix * matB0 = B0->ParallelAssemble(); delete B0;
206 HypreParMatrix * matBhat = Bhat->ParallelAssemble(); delete Bhat;
207 HypreParMatrix * matSinv = Sinv->ParallelAssemble(); delete Sinv;
208 HypreParMatrix * matS0 = S0->ParallelAssemble(); delete S0;
209
210 // 9. Define the block structure of the problem, by creating the offset
211 // variables. Also allocate two BlockVector objects to store the solution
212 // and rhs.
213 enum {x0_var, xhat_var, NVAR};
214
215 int true_s0 = x0_space->TrueVSize();
216 int true_s1 = xhat_space->TrueVSize();
217 int true_s_test = test_space->TrueVSize();
218
219 Array<int> true_offsets(NVAR+1);
220 true_offsets[0] = 0;
221 true_offsets[1] = true_s0;
222 true_offsets[2] = true_s0+true_s1;
223
224 Array<int> true_offsets_test(2);
225 true_offsets_test[0] = 0;
226 true_offsets_test[1] = true_s_test;
227
228 BlockVector x(true_offsets), b(true_offsets);
229 x = 0.0;
230 b = 0.0;
231
232 // 10. Set up the 1x2 block Least Squares DPG operator, B = [B0 Bhat],
233 // the normal equation operator, A = B^t Sinv B, and
234 // the normal equation right-hand-size, b = B^t Sinv F.
235 BlockOperator B(true_offsets_test, true_offsets);
236 B.SetBlock(0, 0, matB0);
237 B.SetBlock(0, 1, matBhat);
238
239 RAPOperator A(B, *matSinv, B);
240
241 HypreParVector *trueF = F->ParallelAssemble();
242 {
243 HypreParVector SinvF(test_space);
244 matSinv->Mult(*trueF, SinvF);
245 B.MultTranspose(SinvF, b);
246 }
247
248 // 11. Set up a block-diagonal preconditioner for the 2x2 normal equation
249 //
250 // [ S0^{-1} 0 ]
251 // [ 0 Shat^{-1} ] Shat = (Bhat^T Sinv Bhat)
252 //
253 // corresponding to the primal (x0) and interfacial (xhat) unknowns.
254 // Since the Shat operator is equivalent to an H(div) matrix reduced to
255 // the interfacial skeleton, we approximate its inverse with one V-cycle
256 // of the ADS preconditioner from the hypre library (in 2D we use AMS for
257 // the rotated H(curl) problem).
258 HypreBoomerAMG *S0inv = new HypreBoomerAMG(*matS0);
259 S0inv->SetPrintLevel(0);
260
261 HypreParMatrix *Shat = RAP(matSinv, matBhat);
262 HypreSolver *Shatinv;
263 if (dim == 2) { Shatinv = new HypreAMS(*Shat, xhat_space); }
264 else { Shatinv = new HypreADS(*Shat, xhat_space); }
265
266 BlockDiagonalPreconditioner P(true_offsets);
267 P.SetDiagonalBlock(0, S0inv);
268 P.SetDiagonalBlock(1, Shatinv);
269
270 // 12. Solve the normal equation system using the PCG iterative solver.
271 // Check the weighted norm of residual for the DPG least square problem.
272 // Wrap the primal variable in a GridFunction for visualization purposes.
273 CGSolver pcg(MPI_COMM_WORLD);
274 pcg.SetOperator(A);
275 pcg.SetPreconditioner(P);
276 pcg.SetRelTol(1e-6);
277 pcg.SetMaxIter(200);
278 pcg.SetPrintLevel(1);
279 pcg.Mult(b, x);
280
281 {
282 HypreParVector LSres(test_space), tmp(test_space);
283 B.Mult(x, LSres);
284 LSres -= *trueF;
285 matSinv->Mult(LSres, tmp);
286 real_t res = sqrt(InnerProduct(LSres, tmp));
287 if (myid == 0)
288 {
289 cout << "\n|| B0*x0 + Bhat*xhat - F ||_{S^-1} = " << res << endl;
290 }
291 }
292
293 x0->Distribute(x.GetBlock(x0_var));
294
295 // 13. Save the refined mesh and the solution in parallel. This output can
296 // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
297 {
298 ostringstream mesh_name, sol_name;
299 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
300 sol_name << "sol." << setfill('0') << setw(6) << myid;
301
302 ofstream mesh_ofs(mesh_name.str().c_str());
303 mesh_ofs.precision(8);
304 pmesh->Print(mesh_ofs);
305
306 ofstream sol_ofs(sol_name.str().c_str());
307 sol_ofs.precision(8);
308 x0->Save(sol_ofs);
309 }
310
311 // 14. Send the solution by socket to a GLVis server.
312 if (visualization)
313 {
314 char vishost[] = "localhost";
315 int visport = 19916;
316 socketstream sol_sock(vishost, visport);
317 sol_sock << "parallel " << num_procs << " " << myid << "\n";
318 sol_sock.precision(8);
319 sol_sock << "solution\n" << *pmesh << *x0 << flush;
320 }
321
322 // 15. Free the used memory.
323 delete trueF;
324 delete Shatinv;
325 delete S0inv;
326 delete Shat;
327 delete matB0;
328 delete matBhat;
329 delete matSinv;
330 delete matS0;
331 delete x0;
332 delete F;
333 delete test_space;
334 delete xhat_space;
335 delete x0_space;
336 delete test_fec;
337 delete xhat_fec;
338 delete x0_fec;
339 delete pmesh;
340
341 return 0;
342}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void EliminateEssentialBC(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate essential boundary DOFs from the system.
A class to handle Block diagonal preconditioners in a matrix-free implementation.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
A class to handle Block systems in a matrix-free implementation.
void SetBlock(int iRow, int iCol, Operator *op, real_t c=1.0)
Add a block op in the block-entry (iblock, jblock).
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
A class to handle Vectors in a block fashion.
Conjugate gradient method.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
A coefficient that is constant across space and time.
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
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The Auxiliary-space Divergence Solver in hypre.
Definition hypre.hpp:1922
The Auxiliary-space Maxwell Solver in hypre.
Definition hypre.hpp:1845
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, real_t alpha=1.0, real_t beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition hypre.cpp:1815
Wrapper for hypre's parallel vector class.
Definition hypre.hpp:206
Abstract class for hypre's solvers and preconditioners.
Definition hypre.hpp:1162
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Integrator that inverts the matrix assembled by another integrator.
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
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
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
int MeshGenerator() const
Get the mesh generator/type.
Definition mesh.hpp:1187
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY.
void Assemble(int skip_zeros=1)
void EliminateEssentialBCFromTrialDofs(const Array< int > &marked_vdofs, const Vector &sol, Vector &rhs)
Eliminate the list of DOFs from the columns of the system.
void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi)
Add a trace face integrator. Assumes ownership of bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
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.
Class for parallel bilinear form.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs, int component=-1) const override
Determine the boundary degrees of freedom.
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition pfespace.hpp:436
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
void Distribute(const Vector *tv)
Class for parallel linear form.
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Class for parallel meshes.
Definition pmesh.hpp:34
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Class for parallel bilinear form using different test and trial FE spaces.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
The operator x -> R*A*P*x constructed through the actions of R^T, A and P.
Definition operator.hpp:817
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:445
Integrator defining a sum of multiple Integrators.
int dim
Definition ex24.cpp:53
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
const int visport
real_t InnerProduct(HypreParVector *x, HypreParVector *y)
Definition hypre.cpp:439
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
float real_t
Definition config.hpp:43
const char vishost[]