1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 //                      -------------------------------
13 //                      Convergence Rates Test (Serial)
14 //                      -------------------------------
15 //
16 // Compile with: make rates
17 //
18 // Sample runs:  rates -m ../../data/inline-segment.mesh -sr 4 -prob 0 -o 1
19 //               rates -m ../../data/inline-quad.mesh -sr 3 -prob 0 -o 2
20 //               rates -m ../../data/inline-quad.mesh -sr 3 -prob 1 -o 2
21 //               rates -m ../../data/inline-quad.mesh -sr 3 -prob 2 -o 2
22 //               rates -m ../../data/inline-tri.mesh -sr 2 -prob 2 -o 3
23 //               rates -m ../../data/star.mesh -sr 2 -prob 1 -o 4
24 //               rates -m ../../data/fichera.mesh -sr 3 -prob 2 -o 1
25 //               rates -m ../../data/inline-wedge.mesh -sr 1 -prob 0 -o 2
26 //               rates -m ../../data/inline-hex.mesh -sr 1 -prob 1 -o 2
27 //               rates -m ../../data/square-disc.mesh -sr 2 -prob 1 -o 1
28 //               rates -m ../../data/star.mesh -sr 2 -prob 3 -o 2
29 //               rates -m ../../data/star.mesh -sr 2 -prob 3 -o 2 -j 0
30 //               rates -m ../../data/inline-hex.mesh -sr 1 -prob 3 -o 1
31 //
32 // Description:  This example code demonstrates the use of MFEM to define and
33 //               solve finite element problem for various discretizations and
34 //               provide convergence rates in serial.
35 //
36 //               prob 0: H1 projection:
37 //                       (grad u, grad v) + (u,v) = (grad u_exact, grad v) + (u_exact, v)
38 //               prob 1: H(curl) projection
39 //                       (curl u, curl v) + (u,v) = (curl u_exact, curl v) + (u_exact, v)
40 //               prob 2: H(div) projection
41 //                       (div  u, div  v) + (u,v) = (div  u_exact, div  v) + (u_exact, v)
42 //               prob 3: DG discretization for the Poisson problem
43 //                       -Delta u = f
44 
45 #include "mfem.hpp"
46 #include <fstream>
47 #include <iostream>
48 using namespace std;
49 using namespace mfem;
50 
51 // Exact solution parameters:
52 double sol_s[3] = { -0.32, 0.15, 0.24 };
53 double sol_k[3] = { 1.21, 1.45, 1.37 };
54 
55 // H1
56 double scalar_u_exact(const Vector &x);
57 double rhs_func(const Vector &x);
58 void gradu_exact(const Vector &x, Vector &gradu);
59 
60 // Vector FE
61 void vector_u_exact(const Vector &x, Vector & vector_u);
62 // H(curl)
63 void curlu_exact(const Vector &x, Vector &curlu);
64 // H(div)
65 double divu_exact(const Vector &x);
66 
67 int dim;
68 int prob=0;
69 
main(int argc,char * argv[])70 int main(int argc, char *argv[])
71 {
72    // 1. Parse command-line options.
73    const char *mesh_file = "../../data/inline-quad.mesh";
74    int order = 1;
75    bool visualization = 1;
76    int sr = 1;
77    int jump_scaling_type = 1;
78    double sigma = -1.0;
79    double kappa = -1.0;
80    OptionsParser args(argc, argv);
81    args.AddOption(&mesh_file, "-m", "--mesh",
82                   "Mesh file to use.");
83    args.AddOption(&order, "-o", "--order",
84                   "Finite element order (polynomial degree)");
85    args.AddOption(&prob, "-prob", "--problem",
86                   "Problem kind: 0: H1, 1: H(curl), 2: H(div), 3: DG ");
87    args.AddOption(&sigma, "-s", "--sigma",
88                   "One of the two DG penalty parameters, typically +1/-1."
89                   " See the documentation of class DGDiffusionIntegrator.");
90    args.AddOption(&kappa, "-k", "--kappa",
91                   "One of the two DG penalty parameters, should be positive."
92                   " Negative values are replaced with (order+1)^2.");
93    args.AddOption(&jump_scaling_type, "-j", "--jump-scaling",
94                   "Scaling of the jump error for DG methods: "
95                   "0: no scaling, 1: 1/h, 2: p^2/h");
96    args.AddOption(&sr, "-sr", "--serial_ref",
97                   "Number of serial refinements.");
98    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
99                   "--no-visualization",
100                   "Enable or disable GLVis visualization.");
101    args.Parse();
102    if (!args.Good())
103    {
104       args.PrintUsage(cout);
105       return 1;
106    }
107    if (prob >3 || prob <0) prob = 0; // default problem = H1
108    if (prob == 3)
109    {
110       if (kappa < 0)
111       {
112          kappa = (order+1)*(order+1);
113       }
114    }
115    args.PrintOptions(cout);
116 
117    // 2. Read the (serial) mesh from the given mesh file.
118    Mesh *mesh = new Mesh(mesh_file, 1, 1);
119    dim = mesh->Dimension();
120 
121    // 3. Refine the serial mesh on all processors to increase the resolution.
122    mesh->UniformRefinement();
123 
124    // 4. Define a finite element space on the parallel mesh.
125    FiniteElementCollection *fec=nullptr;
126    switch (prob)
127    {
128       case 0: fec = new H1_FECollection(order,dim);   break;
129       case 1: fec = new ND_FECollection(order,dim);   break;
130       case 2: fec = new RT_FECollection(order-1,dim); break;
131       case 3: fec = new DG_FECollection(order,dim); break;
132       default: break;
133    }
134    FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
135 
136    // 5. Define the solution vector x as a parallel finite element grid function
137    //    corresponding to fespace.
138    GridFunction x(fespace);
139    x = 0.0;
140 
141    // 6. Set up the linear form b(.) and the bilinear form a(.,.).
142    FunctionCoefficient *f=nullptr;
143    FunctionCoefficient *scalar_u=nullptr;
144    FunctionCoefficient *divu=nullptr;
145    VectorFunctionCoefficient *vector_u=nullptr;
146    VectorFunctionCoefficient *gradu=nullptr;
147    VectorFunctionCoefficient *curlu=nullptr;
148 
149    ConstantCoefficient one(1.0);
150    LinearForm b(fespace);
151    BilinearForm a(fespace);
152 
153    switch (prob)
154    {
155       case 0:
156          //(grad u_ex, grad v) + (u_ex,v)
157          scalar_u = new FunctionCoefficient(scalar_u_exact);
158          gradu = new VectorFunctionCoefficient(dim,gradu_exact);
159          b.AddDomainIntegrator(new DomainLFGradIntegrator(*gradu));
160          b.AddDomainIntegrator(new DomainLFIntegrator(*scalar_u));
161 
162          // (grad u, grad v) + (u,v)
163          a.AddDomainIntegrator(new DiffusionIntegrator(one));
164          a.AddDomainIntegrator(new MassIntegrator(one));
165 
166          break;
167       case 1:
168          //(curl u_ex, curl v) + (u_ex,v)
169          vector_u = new VectorFunctionCoefficient(dim,vector_u_exact);
170          curlu = new VectorFunctionCoefficient((dim==3)?dim:1,curlu_exact);
171          b.AddDomainIntegrator(new VectorFEDomainLFCurlIntegrator(*curlu));
172          b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(*vector_u));
173 
174          // (curl u, curl v) + (u,v)
175          a.AddDomainIntegrator(new CurlCurlIntegrator(one));
176          a.AddDomainIntegrator(new VectorFEMassIntegrator(one));
177          break;
178 
179       case 2:
180          //(div u_ex, div v) + (u_ex,v)
181          vector_u = new VectorFunctionCoefficient(dim,vector_u_exact);
182          divu = new FunctionCoefficient(divu_exact);
183          b.AddDomainIntegrator(new VectorFEDomainLFDivIntegrator(*divu));
184          b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(*vector_u));
185 
186          // (div u, div v) + (u,v)
187          a.AddDomainIntegrator(new DivDivIntegrator(one));
188          a.AddDomainIntegrator(new VectorFEMassIntegrator(one));
189          break;
190 
191       case 3:
192          scalar_u = new FunctionCoefficient(scalar_u_exact);
193          f = new FunctionCoefficient(rhs_func);
194          gradu = new VectorFunctionCoefficient(dim,gradu_exact);
195          b.AddDomainIntegrator(new DomainLFIntegrator(*f));
196          b.AddBdrFaceIntegrator(
197       new DGDirichletLFIntegrator(*scalar_u, one, sigma, kappa));
198          a.AddDomainIntegrator(new DiffusionIntegrator(one));
199          a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
200          a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
201          break;
202 
203       default:
204          break;
205    }
206 
207    // 7. Perform successive refinements, compute the errors and the
208    //    corresponding rates of convergence.
209    ConvergenceStudy rates;
210    for (int l = 0; l <= sr; l++)
211    {
212       b.Assemble();
213       a.Assemble();
214       a.Finalize();
215       const SparseMatrix &A = a.SpMat();
216       GSSmoother M(A);
217       if (prob == 3 && sigma != -1.0)
218       {
219          GMRES(A, M, b, x, 0, 500, 10, 1e-12, 0.0);
220       }
221       else
222       {
223          PCG(A, M, b, x, 0, 500, 1e-12, 0.0);
224       }
225 
226       JumpScaling js(1.0, jump_scaling_type == 2 ? JumpScaling::P_SQUARED_OVER_H
227                         : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
228                         : JumpScaling::CONSTANT);
229 
230       switch (prob)
231       {
232          case 0: rates.AddH1GridFunction(&x,scalar_u,gradu); break;
233          case 1: rates.AddHcurlGridFunction(&x,vector_u,curlu); break;
234          case 2: rates.AddHdivGridFunction(&x,vector_u,divu);  break;
235          case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one,js); break;
236       }
237 
238       if (l==sr) break;
239 
240       mesh->UniformRefinement();
241       fespace->Update();
242       a.Update();
243       b.Update();
244       x.Update();
245    }
246    rates.Print();
247 
248    // 8. Send the solution by socket to a GLVis server.
249    if (visualization)
250    {
251       char vishost[] = "localhost";
252       int  visport   = 19916;
253       socketstream sol_sock(vishost, visport);
254       sol_sock.precision(8);
255       sol_sock << "solution\n" << *mesh << x <<
256                "window_title 'Numerical Pressure (real part)' "
257                << flush;
258    }
259 
260    // 9. Free the used memory.
261    delete f;
262    delete scalar_u;
263    delete divu;
264    delete vector_u;
265    delete gradu;
266    delete curlu;
267    delete fespace;
268    delete fec;
269    delete mesh;
270    return 0;
271 }
272 
rhs_func(const Vector & x)273 double rhs_func(const Vector &x)
274 {
275    double val = 1.0, lap = 0.0;
276    for (int d = 0; d < x.Size(); d++)
277    {
278       const double f = sin(M_PI*(sol_s[d]+sol_k[d]*x(d)));
279       val *= f;
280       lap = lap*f + val*M_PI*M_PI*sol_k[d]*sol_k[d];
281    }
282    return lap;
283 }
284 
scalar_u_exact(const Vector & x)285 double scalar_u_exact(const Vector &x)
286 {
287    double val = 1.0;
288    for (int d = 0; d < x.Size(); d++)
289    {
290       val *= sin(M_PI*(sol_s[d]+sol_k[d]*x(d)));
291    }
292    return val;
293 }
294 
gradu_exact(const Vector & x,Vector & grad)295 void gradu_exact(const Vector &x, Vector &grad)
296 {
297    grad.SetSize(x.Size());
298    double *g = grad.GetData();
299    double val = 1.0;
300    for (int d = 0; d < x.Size(); d++)
301    {
302       const double y = M_PI*(sol_s[d]+sol_k[d]*x(d));
303       const double f = sin(y);
304       for (int j = 0; j < d; j++) { g[j] *= f; }
305       g[d] = val*M_PI*sol_k[d]*cos(y);
306       val *= f;
307    }
308 }
309 
vector_u_exact(const Vector & x,Vector & vector_u)310 void vector_u_exact(const Vector &x, Vector & vector_u)
311 {
312    vector_u.SetSize(x.Size());
313    vector_u=0.0;
314    vector_u[0] = scalar_u_exact(x);
315 }
316 
317 // H(curl)
curlu_exact(const Vector & x,Vector & curlu)318 void curlu_exact(const Vector &x, Vector &curlu)
319 {
320    Vector grad;
321    gradu_exact(x,grad);
322    int n = (x.Size()==3)?3:1;
323    curlu.SetSize(n);
324    if (x.Size()==3)
325    {
326       curlu[0] = 0.0;
327       curlu[1] = grad[2];
328       curlu[2] = -grad[1];
329    }
330    else if (x.Size()==2)
331    {
332       curlu[0] = -grad[1];
333    }
334 }
335 
336 // H(div)
divu_exact(const Vector & x)337 double divu_exact(const Vector &x)
338 {
339    Vector grad;
340    gradu_exact(x,grad);
341 
342    return grad[0];
343 }
344