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