MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
convergence.cpp
Go to the documentation of this file.
1 #include "convergence.hpp"
2 
3 using namespace std;
4 
5 
6 namespace mfem
7 {
8 
9 void ConvergenceStudy::Reset()
10 {
11  counter=0;
12  dcounter=0;
13  fcounter=0;
14  cont_type=-1;
15  print_flag=1;
16  L2Errors.SetSize(0);
17  L2Rates.SetSize(0);
18  DErrors.SetSize(0);
19  DRates.SetSize(0);
20  EnErrors.SetSize(0);
21  EnRates.SetSize(0);
22  DGFaceErrors.SetSize(0);
23  DGFaceRates.SetSize(0);
24  ndofs.SetSize(0);
25 }
26 
27 double ConvergenceStudy::GetNorm(GridFunction *gf, Coefficient *scalar_u,
28  VectorCoefficient *vector_u)
29 {
30  bool norm_set = false;
31  double norm=0.0;
32  int order = gf->FESpace()->GetMaxElementOrder();
33  int order_quad = std::max(2, 2*order+1);
34  const IntegrationRule *irs[Geometry::NumGeom];
35  for (int i=0; i < Geometry::NumGeom; ++i)
36  {
37  irs[i] = &(IntRules.Get(i, order_quad));
38  }
39 #ifdef MFEM_USE_MPI
40  ParGridFunction *pgf = dynamic_cast<ParGridFunction *>(gf);
41  if (pgf)
42  {
43  ParMesh *pmesh = pgf->ParFESpace()->GetParMesh();
44  if (scalar_u)
45  {
46  norm = ComputeGlobalLpNorm(2.0,*scalar_u,*pmesh,irs);
47  }
48  else if (vector_u)
49  {
50  norm = ComputeGlobalLpNorm(2.0,*vector_u,*pmesh,irs);
51  }
52  norm_set = true;
53  }
54 #endif
55  if (!norm_set)
56  {
57  Mesh *mesh = gf->FESpace()->GetMesh();
58  if (scalar_u)
59  {
60  norm = ComputeLpNorm(2.0,*scalar_u,*mesh,irs);
61  }
62  else if (vector_u)
63  {
64  norm = ComputeLpNorm(2.0,*vector_u,*mesh,irs);
65  }
66  }
67  return norm;
68 }
69 
70 void ConvergenceStudy::AddL2Error(GridFunction *gf,
71  Coefficient *scalar_u, VectorCoefficient *vector_u)
72 {
73  int tdofs=0;
74  int dim = gf->FESpace()->GetMesh()->Dimension();
75 
76 #ifdef MFEM_USE_MPI
77  ParGridFunction *pgf = dynamic_cast<ParGridFunction *>(gf);
78  if (pgf)
79  {
80  MPI_Comm comm = pgf->ParFESpace()->GetComm();
81  int rank;
82  MPI_Comm_rank(comm, &rank);
83  print_flag = 0;
84  if (rank==0) { print_flag = 1; }
85  tdofs = pgf->ParFESpace()->GlobalTrueVSize();
86  }
87 #endif
88  if (!tdofs) { tdofs = gf->FESpace()->GetTrueVSize(); }
89  ndofs.Append(tdofs);
90  double L2Err = 1.;
91  if (scalar_u)
92  {
93  L2Err = gf->ComputeL2Error(*scalar_u);
94  CoeffNorm = GetNorm(gf,scalar_u,nullptr);
95  }
96  else if (vector_u)
97  {
98  L2Err = gf->ComputeL2Error(*vector_u);
99  CoeffNorm = GetNorm(gf,nullptr,vector_u);
100  }
101  else
102  {
103  MFEM_ABORT("Exact Solution Coefficient pointer is NULL");
104  }
105  L2Errors.Append(L2Err);
106  // Compute the rate of convergence by:
107  // rate = log (||u - u_h|| / ||u - u_{h/2}||)/(1/dim * log(N_{h/2}/N_{h}))
108  double val=0.;
109  if (counter)
110  {
111  double num = log(L2Errors[counter-1]/L2Err);
112  double den = log((double)ndofs[counter]/ndofs[counter-1]);
113  val = dim * num/den;
114  }
115  L2Rates.Append(val);
116  counter++;
117 }
118 
119 void ConvergenceStudy::AddGf(GridFunction *gf, Coefficient *scalar_u,
120  VectorCoefficient *grad,
121  Coefficient *ell_coeff,
122  JumpScaling jump_scaling)
123 {
124  cont_type = gf->FESpace()->FEColl()->GetContType();
125 
126  MFEM_VERIFY((cont_type == mfem::FiniteElementCollection::CONTINUOUS) ||
127  (cont_type == mfem::FiniteElementCollection::DISCONTINUOUS),
128  "This constructor is intended for H1 or L2 Elements")
129 
130  AddL2Error(gf,scalar_u, nullptr);
131  int dim = gf->FESpace()->GetMesh()->Dimension();
132 
133  if (grad)
134  {
135  double GradErr = gf->ComputeGradError(grad);
136  DErrors.Append(GradErr);
137  double error =
138  sqrt(L2Errors[counter-1]*L2Errors[counter-1]+GradErr*GradErr);
139  EnErrors.Append(error);
140  // Compute the rate of convergence by:
141  // rate = log (||u - u_h|| / ||u - u_{h/2}||)/(1/dim * log(N_{h/2}/N_{h}))
142  double val = 0.;
143  double eval = 0.;
144  if (dcounter)
145  {
146  double num = log(DErrors[dcounter-1]/GradErr);
147  double den = log((double)ndofs[dcounter]/ndofs[dcounter-1]);
148  val = dim * num/den;
149  num = log(EnErrors[dcounter-1]/error);
150  eval = dim * num/den;
151  }
152  DRates.Append(val);
153  EnRates.Append(eval);
154  CoeffDNorm = GetNorm(gf,nullptr,grad);
155  dcounter++;
156  MFEM_VERIFY(counter == dcounter,
157  "Number of added solutions and derivatives do not match")
158  }
159 
160  if (cont_type == mfem::FiniteElementCollection::DISCONTINUOUS && ell_coeff)
161  {
162  double DGErr = gf->ComputeDGFaceJumpError(scalar_u,ell_coeff,jump_scaling);
163  DGFaceErrors.Append(DGErr);
164  // Compute the rate of convergence by:
165  // rate = log (||u - u_h|| / ||u - u_{h/2}||)/(1/dim * log(N_{h/2}/N_{h}))
166  double val = 0.;
167  if (fcounter)
168  {
169  double num = log(DGFaceErrors[fcounter-1]/DGErr);
170  double den = log((double)ndofs[fcounter]/ndofs[fcounter-1]);
171  val = dim * num/den;
172  }
173  DGFaceRates.Append(val);
174  fcounter++;
175  MFEM_VERIFY(fcounter == counter, "Number of added solutions mismatch");
176  }
177 }
178 
179 void ConvergenceStudy::AddGf(GridFunction *gf, VectorCoefficient *vector_u,
180  VectorCoefficient *curl, Coefficient *div)
181 {
182  cont_type = gf->FESpace()->FEColl()->GetContType();
183 
184  AddL2Error(gf,nullptr,vector_u);
185  int dim = gf->FESpace()->GetMesh()->Dimension();
186  double DErr = 0.0;
187  bool derivative = false;
188  if (curl)
189  {
190  DErr = gf->ComputeCurlError(curl);
191  CoeffDNorm = GetNorm(gf,nullptr,curl);
192  derivative = true;
193  }
194  else if (div)
195  {
196  DErr = gf->ComputeDivError(div);
197  // update coefficient norm
198  CoeffDNorm = GetNorm(gf,div,nullptr);
199  derivative = true;
200  }
201  if (derivative)
202  {
203  double error = sqrt(L2Errors[counter-1]*L2Errors[counter-1] + DErr*DErr);
204  DErrors.Append(DErr);
205  EnErrors.Append(error);
206  // Compute the rate of convergence by:
207  // rate = log (||u - u_h|| / ||u - u_{h/2}||)/(1/dim * log(N_{h/2}/N_{h}))
208  double val = 0.;
209  double eval = 0.;
210  if (dcounter)
211  {
212  double num = log(DErrors[dcounter-1]/DErr);
213  double den = log((double)ndofs[dcounter]/ndofs[dcounter-1]);
214  val = dim * num/den;
215  num = log(EnErrors[dcounter-1]/error);
216  eval = dim * num/den;
217  }
218  DRates.Append(val);
219  EnRates.Append(eval);
220  dcounter++;
221  MFEM_VERIFY(counter == dcounter,
222  "Number of added solutions and derivatives do not match")
223  }
224 }
225 
226 void ConvergenceStudy::Print(bool relative, std::ostream &os)
227 {
228  if (print_flag)
229  {
230  std::string title = (relative) ? "Relative " : "Absolute ";
231  os << "\n";
232  os << " -------------------------------------------" << "\n";
233  os << std::setw(21) << title << "L2 Error " << "\n";
234  os << " -------------------------------------------"
235  << "\n";
236  os << std::right<< std::setw(11)<< "DOFs "<< std::setw(13) << "Error ";
237  os << std::setw(15) << "Rate " << "\n";
238  os << " -------------------------------------------"
239  << "\n";
240  os << std::setprecision(4);
241  double d = (relative) ? CoeffNorm : 1.0;
242  for (int i =0; i<counter; i++)
243  {
244  os << std::right << std::setw(10)<< ndofs[i] << std::setw(16)
245  << std::scientific << L2Errors[i]/d << std::setw(13)
246  << std::fixed << L2Rates[i] << "\n";
247  }
248  os << "\n";
249  if (dcounter == counter)
250  {
251  std::string dname;
252  switch (cont_type)
253  {
254  case 0: dname = "Grad"; break;
255  case 1: dname = "Curl"; break;
256  case 2: dname = "Div"; break;
257  case 3: dname = "DG Grad"; break;
258  default: break;
259  }
260  os << " -------------------------------------------" << "\n";
261  os << std::setw(21) << title << dname << " Error " << "\n";
262  os << " -------------------------------------------" << "\n";
263  os << std::right<<std::setw(11)<< "DOFs "<< std::setw(13) << "Error";
264  os << std::setw(15) << "Rate " << "\n";
265  os << " -------------------------------------------"
266  << "\n";
267  os << std::setprecision(4);
268  d = (relative) ? CoeffDNorm : 1.0;
269  for (int i =0; i<dcounter; i++)
270  {
271  os << std::right << std::setw(10)<< ndofs[i] << std::setw(16)
272  << std::scientific << DErrors[i]/d << std::setw(13)
273  << std::fixed << DRates[i] << "\n";
274  }
275  os << "\n";
276  switch (cont_type)
277  {
278  case 0: dname = "H1"; break;
279  case 1: dname = "H(Curl)"; break;
280  case 2: dname = "H(Div)"; break;
281  case 3: dname = "DG H1"; break;
282  default: break;
283  }
284 
285  if (dcounter)
286  {
287  d = (relative) ?
288  sqrt(CoeffNorm*CoeffNorm + CoeffDNorm*CoeffDNorm):1.0;
289 
290  os << " -------------------------------------------" << "\n";
291  os << std::setw(21) << title << dname << " Error " << "\n";
292  os << " -------------------------------------------" << "\n";
293  os << std::right<< std::setw(11)<< "DOFs "<< std::setw(13);
294  os << "Error ";
295  os << std::setw(15) << "Rate " << "\n";
296  os << " -------------------------------------------"
297  << "\n";
298  os << std::setprecision(4);
299  for (int i =0; i<dcounter; i++)
300  {
301  os << std::right << std::setw(10)<< ndofs[i] << std::setw(16)
302  << std::scientific << EnErrors[i]/d << std::setw(13)
303  << std::fixed << EnRates[i] << "\n";
304  }
305  os << "\n";
306  }
307  }
308  if (cont_type == 3 && fcounter)
309  {
310  os << " -------------------------------------------" << "\n";
311  os << " DG Face Jump Error " << "\n";
312  os << " -------------------------------------------"
313  << "\n";
314  os << std::right<< std::setw(11)<< "DOFs "<< std::setw(13);
315  os << "Error ";
316  os << std::setw(15) << "Rate " << "\n";
317  os << " -------------------------------------------"
318  << "\n";
319  os << std::setprecision(4);
320  for (int i =0; i<fcounter; i++)
321  {
322  os << std::right << std::setw(10)<< ndofs[i] << std::setw(16)
323  << std::scientific << DGFaceErrors[i] << std::setw(13)
324  << std::fixed << DGFaceRates[i] << "\n";
325  }
326  os << "\n";
327  }
328  }
329 }
330 
331 } // namespace mfem
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
mfem::GridFunction
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
mfem::IntegrationRules::Get
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
mfem::VectorCoefficient
Base class for vector Coefficients that optionally depend on time and space.
Definition: coefficient.hpp:448
mfem::FiniteElementCollection::DISCONTINUOUS
Field is discontinuous across element interfaces.
Definition: fe_coll.hpp:48
mfem::ComputeLpNorm
double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
Compute the Lp norm of a function f. .
Definition: coefficient.cpp:1416
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
mfem::FiniteElementSpace::GetMaxElementOrder
int GetMaxElementOrder() const
Return the maximum polynomial order.
Definition: fespace.hpp:459
mfem::GridFunction::FESpace
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:668
mfem::Coefficient
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
mfem::ComputeGlobalLpNorm
double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh, const IntegrationRule *irs[])
Compute the global Lp norm of a function f. .
Definition: coefficient.cpp:1459
GetMesh
Mesh * GetMesh(int type)
Definition: ex29.cpp:218
dim
int dim
Definition: ex24.cpp:53
mfem::FiniteElementCollection::CONTINUOUS
Field is continuous across element interfaces.
Definition: fe_coll.hpp:45
mfem::IntRules
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
Generated on Wed Dec 7 2022 10:41:45 for MFEM by   doxygen 1.8.5