MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
diffusion.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// MFEM Ultraweak DPG example for diffusion
13//
14// Compile with: make diffusion
15//
16// sample runs
17// diffusion -m ../../data/star.mesh -o 3 -ref 1 -do 1 -prob 1 -sc
18// diffusion -m ../../data/inline-tri.mesh -o 2 -ref 2 -do 1 -prob 0
19// diffusion -m ../../data/inline-quad.mesh -o 4 -ref 1 -do 2 -prob 0 -sc
20// diffusion -m ../../data/inline-tet.mesh -o 3 -ref 0 -do 1 -prob 1 -sc
21
22// Description:
23// This example code demonstrates the use of MFEM to define and solve
24// the "ultraweak" (UW) DPG formulation for the Poisson problem
25
26// - Δ u = f, in Ω
27// u = u₀, on ∂Ω
28
29// It solves two kinds of problems
30// a) f = 1 and u₀ = 0 (like ex1)
31// b) A manufactured solution problem where u_exact = sin(π * (x + y + z)).
32// This example computes and prints out convergence rates for the L2 error.
33
34// The DPG UW deals with the First Order System
35// ∇ u - σ = 0, in Ω
36// - ∇⋅σ = f, in Ω
37// u = u₀, in ∂Ω
38
39// Ultraweak-DPG is obtained by integration by parts of both equations and the
40// introduction of trace unknowns on the mesh skeleton
41//
42// u ∈ L²(Ω), σ ∈ (L²(Ω))ᵈⁱᵐ
43// û ∈ H^1/2(Γₕ), σ̂ ∈ H^-1/2(Γₕ)
44// -(u , ∇⋅τ) - (σ , τ) + < û, τ⋅n> = 0, ∀ τ ∈ H(div,Ω)
45// (σ , ∇ v) + < σ̂, v > = (f,v) ∀ v ∈ H¹(Ω)
46// û = u₀ on ∂Ω
47
48// Note:
49// û := u and σ̂ := -σ on the mesh skeleton
50//
51// -------------------------------------------------------------
52// | | u | σ | û | σ̂ | RHS |
53// -------------------------------------------------------------
54// | τ | -(u,∇⋅τ) | -(σ,τ) | < û, τ⋅n> | | 0 |
55// | | | | | | |
56// | v | | (σ,∇ v) | | <σ̂,v> | (f,v) |
57
58// where (τ,v) ∈ H(div,Ω) × H^1(Ω)
59
60// Here we use the "space-induced" test norm i.e.,
61//
62// ||(t,v)||²_H(div)×H¹ := ||t||² + ||∇⋅t||² + ||v||² + ||∇v||²
63
64// For more information see https://doi.org/10.1007/978-3-319-01818-8_6
65
66#include "mfem.hpp"
67#include "util/weakform.hpp"
69#include <fstream>
70#include <iostream>
71
72using namespace std;
73using namespace mfem;
74using namespace mfem::common;
75
81
83
84real_t exact_u(const Vector & X);
85void exact_gradu(const Vector & X, Vector &gradu);
87void exact_sigma(const Vector & X, Vector & sigma);
88real_t exact_hatu(const Vector & X);
89void exact_hatsigma(const Vector & X, Vector & hatsigma);
90real_t f_exact(const Vector & X);
91
92int main(int argc, char *argv[])
93{
94 // 1. Parse command-line options.
95 const char *mesh_file = "../../data/inline-quad.mesh";
96 int order = 1;
97 int delta_order = 1;
98 int ref = 0;
99 bool visualization = true;
100 int iprob = 1;
101 bool static_cond = false;
102
103 OptionsParser args(argc, argv);
104 args.AddOption(&mesh_file, "-m", "--mesh",
105 "Mesh file to use.");
106 args.AddOption(&order, "-o", "--order",
107 "Finite element order (polynomial degree).");
108 args.AddOption(&delta_order, "-do", "--delta_order",
109 "Order enrichment for DPG test space.");
110 args.AddOption(&ref, "-ref", "--num_refinements",
111 "Number of uniform refinements");
112 args.AddOption(&iprob, "-prob", "--problem", "Problem case"
113 " 0: manufactured, 1: general");
114 args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
115 "--no-static-condensation", "Enable static condensation.");
116 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
117 "--no-visualization",
118 "Enable or disable GLVis visualization.");
119 args.Parse();
120 if (!args.Good())
121 {
122 args.PrintUsage(cout);
123 return 1;
124 }
125 args.PrintOptions(cout);
126
127 if (iprob > 1) { iprob = 1; }
128 prob = (prob_type)iprob;
129
130 Mesh mesh(mesh_file, 1, 1);
131 int dim = mesh.Dimension();
132 MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example");
133
134 // Define spaces
135 enum TrialSpace
136 {
137 u_space = 0,
138 sigma_space = 1,
139 hatu_space = 2,
140 hatsigma_space = 3
141 };
142 enum TestSpace
143 {
144 tau_space = 0,
145 v_space = 1
146 };
147 // L2 space for u
148 FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
149 FiniteElementSpace *u_fes = new FiniteElementSpace(&mesh,u_fec);
150
151 // Vector L2 space for σ
152 FiniteElementCollection *sigma_fec = new L2_FECollection(order-1,dim);
153 FiniteElementSpace *sigma_fes = new FiniteElementSpace(&mesh,sigma_fec, dim);
154
155 // H^1/2 space for û
156 FiniteElementCollection * hatu_fec = new H1_Trace_FECollection(order,dim);
157 FiniteElementSpace *hatu_fes = new FiniteElementSpace(&mesh,hatu_fec);
158
159 // H^-1/2 space for σ̂
160 FiniteElementCollection * hatsigma_fec = new RT_Trace_FECollection(order-1,dim);
161 FiniteElementSpace *hatsigma_fes = new FiniteElementSpace(&mesh,hatsigma_fec);
162
163 // test space fe collections
164 int test_order = order+delta_order;
165 FiniteElementCollection * tau_fec = new RT_FECollection(test_order-1, dim);
166 FiniteElementCollection * v_fec = new H1_FECollection(test_order, dim);
167
170
171 trial_fes.Append(u_fes);
172 trial_fes.Append(sigma_fes);
173 trial_fes.Append(hatu_fes);
174 trial_fes.Append(hatsigma_fes);
175 test_fec.Append(tau_fec);
176 test_fec.Append(v_fec);
177
178 // Required coefficients for the weak formulation
179 ConstantCoefficient one(1.0);
180 ConstantCoefficient negone(-1.0);
181 FunctionCoefficient f(f_exact); // rhs for the manufactured solution problem
182
183 // Required coefficients for the exact solution case
187
188 // Define the DPG weak formulation
189 DPGWeakForm * a = new DPGWeakForm(trial_fes,test_fec);
190
191 // -(u,∇⋅τ)
192 a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(one),
193 TrialSpace::u_space,TestSpace::tau_space);
194
195 // -(σ,τ)
196 a->AddTrialIntegrator(new TransposeIntegrator(new VectorFEMassIntegrator(
197 negone)), TrialSpace::sigma_space, TestSpace::tau_space);
198
199 // (σ,∇ v)
200 a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(one)),
201 TrialSpace::sigma_space,TestSpace::v_space);
202
203 // <û,τ⋅n>
204 a->AddTrialIntegrator(new NormalTraceIntegrator,
205 TrialSpace::hatu_space,TestSpace::tau_space);
206
207 // -<σ̂,v> (sign is included in σ̂)
208 a->AddTrialIntegrator(new TraceIntegrator,
209 TrialSpace::hatsigma_space, TestSpace::v_space);
210
211 // test integrators (space-induced norm for H(div) × H1)
212 // (∇⋅τ,∇⋅δτ)
213 a->AddTestIntegrator(new DivDivIntegrator(one),
214 TestSpace::tau_space, TestSpace::tau_space);
215 // (τ,δτ)
216 a->AddTestIntegrator(new VectorFEMassIntegrator(one),
217 TestSpace::tau_space, TestSpace::tau_space);
218 // (∇v,∇δv)
219 a->AddTestIntegrator(new DiffusionIntegrator(one),
220 TestSpace::v_space, TestSpace::v_space);
221 // (v,δv)
222 a->AddTestIntegrator(new MassIntegrator(one),
223 TestSpace::v_space, TestSpace::v_space);
224
225 // RHS
226 if (prob == prob_type::manufactured)
227 {
228 a->AddDomainLFIntegrator(new DomainLFIntegrator(f),TestSpace::v_space);
229 }
230 else
231 {
232 a->AddDomainLFIntegrator(new DomainLFIntegrator(one),TestSpace::v_space);
233 }
234
235 // GridFunction for Dirichlet bdr data
236 GridFunction hatu_gf;
237
238 // Visualization streams
239 socketstream u_out;
240 socketstream sigma_out;
241
242 if (prob == prob_type::manufactured)
243 {
244 std::cout << "\n Ref |"
245 << " Dofs |"
246 << " L2 Error |"
247 << " Rate |"
248 << " PCG it |" << endl;
249 std::cout << std::string(50,'-')
250 << endl;
251 }
252
253 real_t err0 = 0.;
254 int dof0=0.;
255 if (static_cond) { a->EnableStaticCondensation(); }
256 for (int it = 0; it<=ref; it++)
257 {
258 a->Assemble();
259
260 Array<int> ess_tdof_list;
261 Array<int> ess_bdr;
262 if (mesh.bdr_attributes.Size())
263 {
264 ess_bdr.SetSize(mesh.bdr_attributes.Max());
265 ess_bdr = 1;
266 hatu_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
267 }
268
269 // shift the ess_tdofs
270 for (int i = 0; i < ess_tdof_list.Size(); i++)
271 {
272 ess_tdof_list[i] += u_fes->GetTrueVSize() + sigma_fes->GetTrueVSize();
273 }
274
275 Array<int> offsets(5);
276 offsets[0] = 0;
277 offsets[1] = u_fes->GetVSize();
278 offsets[2] = sigma_fes->GetVSize();
279 offsets[3] = hatu_fes->GetVSize();
280 offsets[4] = hatsigma_fes->GetVSize();
281 offsets.PartialSum();
282
283 BlockVector x(offsets);
284 x = 0.0;
285 if (prob == prob_type::manufactured)
286 {
287 hatu_gf.MakeRef(hatu_fes,x.GetBlock(2),0);
288 hatu_gf.ProjectBdrCoefficient(hatuex,ess_bdr);
289 }
290
291 OperatorPtr Ah;
292 Vector X,B;
293 a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
294
295 BlockMatrix * A = Ah.As<BlockMatrix>();
296
298 M.owns_blocks = 1;
299 for (int i=0; i<A->NumRowBlocks(); i++)
300 {
301 M.SetDiagonalBlock(i,new GSSmoother(A->GetBlock(i,i)));
302 }
303
304 CGSolver cg;
305 cg.SetRelTol(1e-10);
306 cg.SetMaxIter(2000);
307 cg.SetPrintLevel(prob== prob_type::general ? 3 : 0);
308 cg.SetPreconditioner(M);
309 cg.SetOperator(*A);
310 cg.Mult(B, X);
311
312 a->RecoverFEMSolution(X,x);
313
314 GridFunction u_gf, sigma_gf;
315 u_gf.MakeRef(u_fes,x.GetBlock(0),0);
316 sigma_gf.MakeRef(sigma_fes,x.GetBlock(1),0);
317
318 if (prob == prob_type::manufactured)
319 {
320 int l2dofs = u_fes->GetVSize() + sigma_fes->GetVSize();
321 real_t u_err = u_gf.ComputeL2Error(uex);
322 real_t sigma_err = sigma_gf.ComputeL2Error(sigmaex);
323 real_t L2Error = sqrt(u_err*u_err + sigma_err*sigma_err);
324 real_t rate_err = (it) ? dim*log(err0/L2Error)/log((real_t)dof0/l2dofs) : 0.0;
325 err0 = L2Error;
326 dof0 = l2dofs;
327
328 std::ios oldState(nullptr);
329 oldState.copyfmt(std::cout);
330 std::cout << std::right << std::setw(5) << it << " | "
331 << std::setw(10) << dof0 << " | "
332 << std::setprecision(3)
333 << std::setw(10) << std::scientific << err0 << " | "
334 << std::setprecision(2)
335 << std::setw(6) << std::fixed << rate_err << " | "
336 << std::setw(6) << std::fixed << cg.GetNumIterations() << " | "
337 << std::endl;
338 std::cout.copyfmt(oldState);
339 }
340
341 if (visualization)
342 {
343 const char * keys = (it == 0 && dim == 2) ? "jRcm\n" : nullptr;
344 char vishost[] = "localhost";
345 int visport = 19916;
346 VisualizeField(u_out,vishost, visport, u_gf,
347 "Numerical u", 0,0, 500, 500, keys);
348 VisualizeField(sigma_out,vishost, visport, sigma_gf,
349 "Numerical flux", 500,0,500, 500, keys);
350 }
351
352 if (it == ref) { break; }
353
354 mesh.UniformRefinement();
355 for (int i =0; i<trial_fes.Size(); i++)
356 {
357 trial_fes[i]->Update(false);
358 }
359 a->Update();
360 }
361
362 delete a;
363 delete tau_fec;
364 delete v_fec;
365 delete hatsigma_fes;
366 delete hatsigma_fec;
367 delete hatu_fes;
368 delete hatu_fec;
369 delete sigma_fec;
370 delete sigma_fes;
371 delete u_fec;
372 delete u_fes;
373
374 return 0;
375}
376
378{
379 real_t alpha = M_PI * (X.Sum());
380 return sin(alpha);
381}
382
383void exact_gradu(const Vector & X, Vector & du)
384{
385 du.SetSize(X.Size());
386 real_t alpha = M_PI * (X.Sum());
387 du.SetSize(X.Size());
388 for (int i = 0; i<du.Size(); i++)
389 {
390 du[i] = M_PI * cos(alpha);
391 }
392}
393
395{
396 real_t alpha = M_PI * (X.Sum());
397 real_t u = sin(alpha);
398 return - M_PI*M_PI * u * X.Size();
399}
400
401void exact_sigma(const Vector & X, Vector & sigma)
402{
403 // σ = ∇ u
405}
406
408{
409 return exact_u(X);
410}
411
412void exact_hatsigma(const Vector & X, Vector & hatsigma)
413{
414 exact_sigma(X,hatsigma);
415 hatsigma *= -1.;
416}
417
419{
420 return -exact_laplacian_u(X);
421}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
A class to handle Block diagonal preconditioners in a matrix-free implementation.
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
Array< int > & RowOffsets()
Return the row offsets for block starts.
int NumRowBlocks() const
Return the number of row blocks.
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
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 representing the DPG weak formulation. Given the variational formulation a(u,...
Definition weakform.hpp:34
for Raviart-Thomas elements
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:588
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
A general function coefficient.
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void MakeRef(FiniteElementSpace *f, real_t *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition gridfunc.cpp:203
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
void ProjectBdrCoefficient(Coefficient &coeff, const Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Definition gridfunc.hpp:469
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:322
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition solvers.hpp:260
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
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 Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
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.
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:445
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1283
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
real_t sigma(const Vector &x)
Definition maxwell.cpp:102
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
prob_type
Definition ex25.cpp:149
int main()
real_t a
Definition lissajous.cpp:41
void exact_gradu(const Vector &X, Vector &gradu)
real_t exact_hatu(const Vector &X)
real_t exact_laplacian_u(const Vector &X)
real_t exact_u(const Vector &X)
void exact_sigma(const Vector &X, Vector &sigma)
void exact_hatsigma(const Vector &X, Vector &hatsigma)
prob_type prob
Definition diffusion.cpp:82
prob_type
Definition diffusion.cpp:77
@ general
Definition diffusion.cpp:79
@ manufactured
Definition diffusion.cpp:78
real_t f_exact(const Vector &X)
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]