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 #ifndef MFEM_CONVERGENCE
13 #define MFEM_CONVERGENCE
14 
15 #include "../linalg/linalg.hpp"
16 #include "gridfunc.hpp"
17 #ifdef MFEM_USE_MPI
18 #include "pgridfunc.hpp"
19 #endif
20 
21 namespace mfem
22 {
23 
24 /** @brief Class to compute error and convergence rates.
25     It supports H1, H(curl) (ND elements), H(div) (RT elements) and L2 (DG).
26 
27     For "smooth enough" solutions the Galerkin error measured in the appropriate
28     norm satisfies || u - u_h || ~ h^k
29 
30     Here, k is called the asymptotic rate of convergence
31 
32     For successive uniform h-refinements the rate can be estimated by
33     k = log(||u - u_h|| / ||u - u_{h/2}||)/log(2)
34 */
35 class ConvergenceStudy
36 {
37 private:
38    // counters for solutions/derivatives
39    int counter=0;
40    int dcounter=0;
41    int fcounter=0;
42 
43    // space continuity type
44    int cont_type=-1;
45 
46    // printing flag for helpful for MPI calls
47    int print_flag=1;
48 
49    // exact solution and derivatives
50    double CoeffNorm;
51    double CoeffDNorm;
52 
53    // Arrays to store error/rates
54    Array<double> L2Errors, DGFaceErrors, DErrors, EnErrors;
55    Array<double> L2Rates, DGFaceRates, DRates, EnRates;
56    Array<int> ndofs;
57 
58    void AddL2Error(GridFunction *gf, Coefficient *scalar_u,
59                    VectorCoefficient *vector_u);
60    void AddGf(GridFunction *gf, Coefficient *scalar_u,
61               VectorCoefficient *grad=nullptr,
62               Coefficient *ell_coeff=nullptr,
63               JumpScaling jump_scaling = {1.0, JumpScaling::ONE_OVER_H});
64    void AddGf(GridFunction *gf, VectorCoefficient *vector_u,
65               VectorCoefficient *curl, Coefficient *div);
66    // returns the L2-norm of scalar_u or vector_u
67    double GetNorm(GridFunction *gf, Coefficient *scalar_u,
68                   VectorCoefficient *vector_u);
69 
70 public:
71 
72    /// Clear any internal data
73    void Reset();
74 
75    /// Add L2 GridFunction, the exact solution and possibly its gradient and/or
76    /// DG face jumps parameters
AddL2GridFunction(GridFunction * gf,Coefficient * scalar_u,VectorCoefficient * grad=nullptr,Coefficient * ell_coeff=nullptr,JumpScaling jump_scaling={1.0, JumpScaling::ONE_OVER_H})77    void AddL2GridFunction(GridFunction *gf, Coefficient *scalar_u,
78                           VectorCoefficient *grad=nullptr,
79                           Coefficient *ell_coeff=nullptr,
80                           JumpScaling jump_scaling = {1.0, JumpScaling::ONE_OVER_H})
81    {
82       AddGf(gf, scalar_u, grad, ell_coeff, jump_scaling);
83    }
84 
85    /// Add H1 GridFunction, the exact solution and possibly its gradient
AddH1GridFunction(GridFunction * gf,Coefficient * scalar_u,VectorCoefficient * grad=nullptr)86    void AddH1GridFunction(GridFunction *gf, Coefficient *scalar_u,
87                           VectorCoefficient *grad=nullptr)
88    {
89       AddGf(gf, scalar_u, grad);
90    }
91 
92    /// Add H(curl) GridFunction, the exact solution and possibly its curl
AddHcurlGridFunction(GridFunction * gf,VectorCoefficient * vector_u,VectorCoefficient * curl=nullptr)93    void AddHcurlGridFunction(GridFunction *gf, VectorCoefficient *vector_u,
94                              VectorCoefficient *curl=nullptr)
95    {
96       AddGf(gf, vector_u, curl, nullptr);
97    }
98 
99    /// Add H(div) GridFunction, the exact solution and possibly its div
AddHdivGridFunction(GridFunction * gf,VectorCoefficient * vector_u,Coefficient * div=nullptr)100    void AddHdivGridFunction(GridFunction *gf, VectorCoefficient *vector_u,
101                             Coefficient *div=nullptr)
102    {
103       AddGf(gf,vector_u, nullptr, div);
104    }
105 
106    /// Get the L2 error at step n
GetL2Error(int n)107    double GetL2Error(int n)
108    {
109       MFEM_VERIFY( n <= counter,"Step out of bounds")
110       return L2Errors[n];
111    }
112 
113    /// Get all L2 errors
GetL2Errors(Array<double> & L2Errors_)114    void GetL2Errors(Array<double> & L2Errors_)
115    {
116       L2Errors_ = L2Errors;
117    }
118 
119    /// Get the Grad/Curl/Div error at step n
GetDError(int n)120    double GetDError(int n)
121    {
122       MFEM_VERIFY(n <= dcounter,"Step out of bounds")
123       return DErrors[n];
124    }
125 
126    /// Get all Grad/Curl/Div errors
GetDErrors(Array<double> & DErrors_)127    void GetDErrors(Array<double> & DErrors_)
128    {
129       DErrors_ = DErrors;
130    }
131 
132    /// Get the DGFaceJumps error at step n
GetDGFaceJumpsError(int n)133    double GetDGFaceJumpsError(int n)
134    {
135       MFEM_VERIFY(n<= fcounter,"Step out of bounds")
136       return DGFaceErrors[n];
137    }
138 
139    /// Get all DGFaceJumps errors
GetDGFaceJumpsErrors(Array<double> & DGFaceErrors_)140    void GetDGFaceJumpsErrors(Array<double> & DGFaceErrors_)
141    {
142       DGFaceErrors_ = DGFaceErrors;
143    }
144 
145    /// Print rates and errors
146    void Print(bool relative = false, std::ostream &out = mfem::out);
147 };
148 
149 } // namespace mfem
150 
151 #endif // MFEM_CONVERGENCE
152