MFEM  v4.6.0
Finite element discretization library
pconvection-diffusion.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 parallel example for convection-diffusion
13 //
14 // Compile with: make pconvection-diffusion
15 //
16 // sample runs
17 // mpirun -np 4 pconvection-diffusion -o 2 -ref 3 -prob 0 -eps 1e-1 -beta '4 2' -theta 0.0
18 // mpirun -np 4 pconvection-diffusion -o 3 -ref 3 -prob 0 -eps 1e-2 -beta '2 3' -theta 0.0
19 // mpirun -np 4 pconvection-diffusion -m ../../data/inline-hex.mesh -o 2 -ref 1 -prob 0 -sc -eps 1e-1 -theta 0.0
20 
21 // AMR runs
22 // mpirun -np 4 pconvection-diffusion -o 3 -ref 15 -prob 1 -eps 1e-3 -beta '1 0' -theta 0.7 -sc
23 // mpirun -np 4 pconvection-diffusion -o 3 -ref 20 -prob 2 -eps 5e-3 -theta 0.7 -sc
24 // mpirun -np 4 pconvection-diffusion -o 2 -ref 20 -prob 3 -eps 1e-2 -beta '1 2' -theta 0.7 -sc
25 
26 // Description:
27 // This example code demonstrates the use of MFEM to define and solve a parallel
28 // "ultraweak" (UW) DPG formulation for the convection-diffusion problem
29 
30 // - εΔu + ∇⋅(βu) = f, in Ω
31 // u = u₀ , on ∂Ω
32 
33 // It solves the following kinds of problems
34 // (a) A manufactured solution where u_exact = sin(π * (x + y + z)).
35 // (b) The 2D Erickson-Johnson problem
36 // (c) Internal layer problem
37 // (d) Boundary layer problem
38 
39 // The DPG UW deals with the First Order System
40 // - ∇⋅σ + ∇⋅(βu) = f, in Ω
41 // 1/ε σ - ∇u = 0, in Ω
42 // u = u₀ , on ∂Ω
43 
44 // Ultraweak-DPG is obtained by integration by parts of both equations and the
45 // introduction of trace unknowns on the mesh skeleton
46 //
47 // u ∈ L²(Ω), σ ∈ (L²(Ω))ᵈⁱᵐ
48 // û ∈ H^1/2, f̂ ∈ H^-1/2
49 // -(βu , ∇v) + (σ , ∇v) + < f̂ , v > = (f,v), ∀ v ∈ H¹(Ω)
50 // (u , ∇⋅τ) + 1/ε (σ , τ) + < û , τ⋅n > = 0, ∀ τ ∈ H(div,Ω)
51 // û = u₀ on ∂Ω
52 
53 // Note:
54 // f̂ := βu - σ, û := -u on the mesh skeleton
55 
56 // -------------------------------------------------------------
57 // | | u | σ | û | f̂ | RHS |
58 // -------------------------------------------------------------
59 // | v |-(βu , ∇v) | (σ , ∇v) | | < f̂ ,v > | (f,v) |
60 // | | | | | | |
61 // | τ | (u ,∇⋅τ) | 1/ε(σ , τ)| <û,τ⋅n> | | 0 |
62 
63 // where (v,τ) ∈ H¹(Ωₕ) × H(div,Ωₕ)
64 
65 // For more information see https://doi.org/10.1016/j.camwa.2013.06.010
66 
67 #include "mfem.hpp"
68 #include "util/pweakform.hpp"
69 #include "../common/mfem-common.hpp"
70 #include <fstream>
71 #include <iostream>
72 
73 using namespace std;
74 using namespace mfem;
75 using namespace mfem::common;
76 
78 {
80  EJ, // see https://doi.org/10.1016/j.camwa.2013.06.010
81  curved_streamlines, // see https://doi.org/10.1515/cmam-2018-0207
82  bdr_layer // see https://doi.org/10.1002/num.20640
83 };
84 
85 static const char *enum_str[] =
86 {
87  "sinusoidal",
88  "EJ",
89  "curved_streamlines",
90  "bdr_layer"
91 };
92 
95 double epsilon;
96 
97 double exact_u(const Vector & X);
98 void exact_gradu(const Vector & X, Vector & du);
99 double exact_laplacian_u(const Vector & X);
100 double exact_u(const Vector & X);
101 void exact_sigma(const Vector & X, Vector & sigma);
102 double exact_hatu(const Vector & X);
103 void exact_hatf(const Vector & X, Vector & hatf);
104 double f_exact(const Vector & X);
105 double bdr_data(const Vector &X);
106 void beta_function(const Vector & X, Vector & beta_val);
108 
109 int main(int argc, char *argv[])
110 {
111  Mpi::Init();
112  int myid = Mpi::WorldRank();
113  Hypre::Init();
114 
115  // 1. Parse command-line options.
116  const char *mesh_file = "../../data/inline-quad.mesh";
117  int order = 1;
118  int delta_order = 1;
119  int ref = 1;
120  int iprob = 0;
121  double theta = 0.7;
122  bool static_cond = false;
123  epsilon = 1e0;
124 
125  bool visualization = true;
126  bool paraview = false;
127 
128  OptionsParser args(argc, argv);
129  args.AddOption(&mesh_file, "-m", "--mesh",
130  "Mesh file to use.");
131  args.AddOption(&order, "-o", "--order",
132  "Finite element order (polynomial degree).");
133  args.AddOption(&delta_order, "-do", "--delta-order",
134  "Order enrichment for DPG test space.");
135  args.AddOption(&epsilon, "-eps", "--epsilon",
136  "Epsilon coefficient");
137  args.AddOption(&ref, "-ref", "--num-refinements",
138  "Number of uniform refinements");
139  args.AddOption(&theta, "-theta", "--theta",
140  "Theta parameter for AMR");
141  args.AddOption(&iprob, "-prob", "--problem", "Problem case"
142  " 0: lshape, 1: General");
143  args.AddOption(&beta, "-beta", "--beta",
144  "Vector Coefficient beta");
145  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
146  "--no-static-condensation", "Enable static condensation.");
147  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
148  "--no-visualization",
149  "Enable or disable GLVis visualization.");
150  args.AddOption(&paraview, "-paraview", "--paraview", "-no-paraview",
151  "--no-paraview",
152  "Enable or disable ParaView visualization.");
153  args.Parse();
154  if (!args.Good())
155  {
156  if (myid == 0)
157  {
158  args.PrintUsage(cout);
159  }
160  return 1;
161  }
162 
163  if (iprob > 3) { iprob = 3; }
164  prob = (prob_type)iprob;
165 
168  {
169  mesh_file = "../../data/inline-quad.mesh";
170  }
171 
172  Mesh mesh(mesh_file, 1, 1);
173  int dim = mesh.Dimension();
174  MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example");
175 
176  bool exact_known = true;
177  switch (prob)
178  {
179  case sinusoidal:
180  case EJ:
181  {
182  if (beta.Size() == 0)
183  {
184  beta.SetSize(dim);
185  beta = 0.0;
186  beta[0] = 1.;
187  }
188  break;
189  }
190  case bdr_layer:
191  {
192  beta.SetSize(dim);
193  beta[0] = 1.;
194  beta[1] = 2.;
195  exact_known = false;
196  }
197  break;
198  default:
199  // do nothing; beta is defined as a FunctionCoefficient
200  break;
201  }
202 
203  if (myid == 0)
204  {
205  args.PrintOptions(cout);
206  }
207 
208  mesh.EnsureNCMesh(true);
209 
210  ParMesh pmesh(MPI_COMM_WORLD, mesh);
211  mesh.Clear();
212 
213  // Define spaces
214  enum TrialSpace
215  {
216  u_space = 0,
217  sigma_space = 1,
218  hatu_space = 2,
219  hatf_space = 3
220  };
221  enum TestSpace
222  {
223  v_space = 0,
224  tau_space = 1
225  };
226  // L2 space for u
227  FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
228  ParFiniteElementSpace *u_fes = new ParFiniteElementSpace(&pmesh,u_fec);
229 
230  // Vector L2 space for σ
231  FiniteElementCollection *sigma_fec = new L2_FECollection(order-1,dim);
232  ParFiniteElementSpace *sigma_fes = new ParFiniteElementSpace(&pmesh,sigma_fec,
233  dim);
234 
235  // H^1/2 space for û
236  FiniteElementCollection * hatu_fec = new H1_Trace_FECollection(order,dim);
237  ParFiniteElementSpace *hatu_fes = new ParFiniteElementSpace(&pmesh,hatu_fec);
238 
239  // H^-1/2 space for σ̂
240  FiniteElementCollection * hatf_fec = new RT_Trace_FECollection(order-1,dim);
241  ParFiniteElementSpace *hatf_fes = new ParFiniteElementSpace(&pmesh,hatf_fec);
242 
243  // testspace fe collections
244  int test_order = order+delta_order;
245  FiniteElementCollection * v_fec = new H1_FECollection(test_order, dim);
246  FiniteElementCollection * tau_fec = new RT_FECollection(test_order-1, dim);
247 
248  // Coefficients
249  ConstantCoefficient one(1.0);
250  ConstantCoefficient negone(-1.0);
252  ConstantCoefficient eps1(1./epsilon);
253  ConstantCoefficient negeps1(-1./epsilon);
255  ConstantCoefficient negeps(-epsilon);
256 
258  ScalarVectorProductCoefficient negbetacoeff(-1.0,betacoeff);
259  OuterProductCoefficient bbtcoeff(betacoeff,betacoeff);
260 
261  // Normal equation weak formulation
264 
265  trial_fes.Append(u_fes);
266  trial_fes.Append(sigma_fes);
267  trial_fes.Append(hatu_fes);
268  trial_fes.Append(hatf_fes);
269  test_fec.Append(v_fec);
270  test_fec.Append(tau_fec);
271 
272  ParDPGWeakForm * a = new ParDPGWeakForm(trial_fes,test_fec);
273  a->StoreMatrices(true);
274 
275  //-(βu , ∇v)
276  a->AddTrialIntegrator(new MixedScalarWeakDivergenceIntegrator(betacoeff),
277  TrialSpace::u_space, TestSpace::v_space);
278 
279  // (σ,∇ v)
280  a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(one)),
281  TrialSpace::sigma_space, TestSpace::v_space);
282 
283  // (u ,∇⋅τ)
284  a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(negone),
285  TrialSpace::u_space, TestSpace::tau_space);
286 
287  // 1/ε (σ,τ)
288  a->AddTrialIntegrator(new TransposeIntegrator(new VectorFEMassIntegrator(eps1)),
289  TrialSpace::sigma_space, TestSpace::tau_space);
290 
291  // <û,τ⋅n>
292  a->AddTrialIntegrator(new NormalTraceIntegrator,
293  TrialSpace::hatu_space, TestSpace::tau_space);
294 
295  // <f̂ ,v>
296  a->AddTrialIntegrator(new TraceIntegrator,
297  TrialSpace::hatf_space, TestSpace::v_space);
298 
299 
300  FiniteElementCollection *coeff_fec = new L2_FECollection(0,dim);
301  ParFiniteElementSpace *coeff_fes = new ParFiniteElementSpace(&pmesh,coeff_fec);
302  ParGridFunction c1_gf, c2_gf;
303  GridFunctionCoefficient c1_coeff(&c1_gf);
304  GridFunctionCoefficient c2_coeff(&c2_gf);
305 
306  c1_gf.SetSpace(coeff_fes);
307  c2_gf.SetSpace(coeff_fes);
308  setup_test_norm_coeffs(c1_gf,c2_gf);
309 
310  // c1 (v,δv)
311  a->AddTestIntegrator(new MassIntegrator(c1_coeff),
312  TestSpace::v_space, TestSpace::v_space);
313  // ε (∇v,∇δv)
314  a->AddTestIntegrator(new DiffusionIntegrator(eps),
315  TestSpace::v_space, TestSpace::v_space);
316  // (β⋅∇v, β⋅∇δv)
317  a->AddTestIntegrator(new DiffusionIntegrator(bbtcoeff),
318  TestSpace::v_space, TestSpace::v_space);
319  // c2 (τ,δτ)
320  a->AddTestIntegrator(new VectorFEMassIntegrator(c2_coeff),
321  TestSpace::tau_space, TestSpace::tau_space);
322  // (∇⋅τ,∇⋅δτ)
323  a->AddTestIntegrator(new DivDivIntegrator(one),
324  TestSpace::tau_space, TestSpace::tau_space);
325 
327  if (prob == prob_type::sinusoidal ||
329  {
330  a->AddDomainLFIntegrator(new DomainLFIntegrator(f),TestSpace::v_space);
331  }
332 
335  Array<int> elements_to_refine;
338 
339  ParGridFunction hatu_gf;
340  ParGridFunction hatf_gf;
341 
342  socketstream u_out;
343  socketstream sigma_out;
344 
345  double res0 = 0.;
346  double err0 = 0.;
347  int dof0 = 0; // init to suppress gcc warning
348  if (myid == 0)
349  {
350  std::cout << " Ref |"
351  << " Dofs |" ;
352  if (exact_known)
353  {
354  mfem::out << " L2 Error |"
355  << " Rate |";
356  }
357  std::cout << " Residual |"
358  << " Rate |"
359  << " CG it |" << endl;
360  std::cout << std::string((exact_known) ? 72 : 50,'-')
361  << endl;
362  }
363 
364  if (static_cond) { a->EnableStaticCondensation(); }
365 
366  ParGridFunction u_gf(u_fes); u_gf = 0.0;
367  ParGridFunction sigma_gf(sigma_fes); sigma_gf = 0.0;
368 
369  ParaViewDataCollection * paraview_dc = nullptr;
370 
371  if (paraview)
372  {
373  paraview_dc = new ParaViewDataCollection(enum_str[prob], &pmesh);
374  paraview_dc->SetPrefixPath("ParaView/Convection-Diffusion");
375  paraview_dc->SetLevelsOfDetail(order);
376  paraview_dc->SetCycle(0);
377  paraview_dc->SetDataFormat(VTKFormat::BINARY);
378  paraview_dc->SetHighOrderOutput(true);
379  paraview_dc->SetTime(0.0); // set the time
380  paraview_dc->RegisterField("u",&u_gf);
381  paraview_dc->RegisterField("sigma",&sigma_gf);
382  }
383 
384  for (int it = 0; it<=ref; it++)
385  {
386  a->Assemble();
387 
388  Array<int> ess_tdof_list_uhat;
389  Array<int> ess_tdof_list_fhat;
390  Array<int> ess_bdr_uhat;
391  Array<int> ess_bdr_fhat;
392  if (pmesh.bdr_attributes.Size())
393  {
394  ess_bdr_uhat.SetSize(pmesh.bdr_attributes.Max());
395  ess_bdr_fhat.SetSize(pmesh.bdr_attributes.Max());
396 
397  if (prob == prob_type::EJ)
398  {
399  ess_bdr_uhat = 0;
400  ess_bdr_fhat = 1;
401  ess_bdr_uhat[1] = 1;
402  ess_bdr_fhat[1] = 0;
403  }
404  else
405  {
406  ess_bdr_uhat = 1;
407  ess_bdr_fhat = 0;
408  }
409 
410  hatu_fes->GetEssentialTrueDofs(ess_bdr_uhat, ess_tdof_list_uhat);
411  hatf_fes->GetEssentialTrueDofs(ess_bdr_fhat, ess_tdof_list_fhat);
412  }
413 
414  // shift the ess_tdofs
415  int n = ess_tdof_list_uhat.Size();
416  int m = ess_tdof_list_fhat.Size();
417  Array<int> ess_tdof_list(n+m);
418  for (int j = 0; j < n; j++)
419  {
420  ess_tdof_list[j] = ess_tdof_list_uhat[j]
421  + u_fes->GetTrueVSize()
422  + sigma_fes->GetTrueVSize();
423  }
424  for (int j = 0; j < m; j++)
425  {
426  ess_tdof_list[j+n] = ess_tdof_list_fhat[j]
427  + u_fes->GetTrueVSize()
428  + sigma_fes->GetTrueVSize()
429  + hatu_fes->GetTrueVSize();
430  }
431 
432  Array<int> offsets(5);
433  offsets[0] = 0;
434  offsets[1] = u_fes->GetVSize();
435  offsets[2] = sigma_fes->GetVSize();
436  offsets[3] = hatu_fes->GetVSize();
437  offsets[4] = hatf_fes->GetVSize();
438  offsets.PartialSum();
439  BlockVector x(offsets);
440  x = 0.0;
441  hatu_gf.MakeRef(hatu_fes,x.GetBlock(2),0);
443  hatu_gf.ProjectBdrCoefficient(bdr_cf,ess_bdr_uhat);
444 
445  hatf_gf.MakeRef(hatf_fes,x.GetBlock(3),0);
446  hatf_gf.ProjectBdrCoefficientNormal(hatfex,ess_bdr_fhat);
447 
448  OperatorPtr Ah;
449  Vector X,B;
450  a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
451 
452  BlockOperator * A = Ah.As<BlockOperator>();
453 
455  M.owns_blocks = 1;
456  int skip = 0;
457  if (!static_cond)
458  {
459  HypreBoomerAMG * amg0 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(0,0));
460  HypreBoomerAMG * amg1 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(1,1));
461  amg0->SetPrintLevel(0);
462  amg1->SetPrintLevel(0);
463  M.SetDiagonalBlock(0,amg0);
464  M.SetDiagonalBlock(1,amg1);
465  skip = 2;
466  }
467  HypreBoomerAMG * amg2 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(skip,
468  skip));
469  amg2->SetPrintLevel(0);
470  M.SetDiagonalBlock(skip,amg2);
471 
472  HypreSolver * prec;
473  if (dim == 2)
474  {
475  // AMS preconditioner for 2D H(div) (trace) space
476  prec = new HypreAMS((HypreParMatrix &)A->GetBlock(skip+1,skip+1), hatf_fes);
477  }
478  else
479  {
480  // ADS preconditioner for 3D H(div) (trace) space
481  prec = new HypreADS((HypreParMatrix &)A->GetBlock(skip+1,skip+1), hatf_fes);
482  }
483  M.SetDiagonalBlock(skip+1,prec);
484 
485  CGSolver cg(MPI_COMM_WORLD);
486  cg.SetRelTol(1e-12);
487  cg.SetMaxIter(2000);
488  cg.SetPrintLevel(0);
489  cg.SetPreconditioner(M);
490  cg.SetOperator(*A);
491  cg.Mult(B, X);
492  int num_iter = cg.GetNumIterations();
493 
494  a->RecoverFEMSolution(X,x);
495  Vector & residuals = a->ComputeResidual(x);
496 
497  double residual = residuals.Norml2();
498  double maxresidual = residuals.Max();
499 
500  double gresidual = residual * residual;
501 
502  MPI_Allreduce(MPI_IN_PLACE,&maxresidual,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
503  MPI_Allreduce(MPI_IN_PLACE,&gresidual,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
504 
505  gresidual = sqrt(gresidual);
506 
507  elements_to_refine.SetSize(0);
508  for (int iel = 0; iel<pmesh.GetNE(); iel++)
509  {
510  if (residuals[iel] > theta * maxresidual)
511  {
512  elements_to_refine.Append(iel);
513  }
514  }
515 
516  u_gf.MakeRef(u_fes,x.GetBlock(0),0);
517  sigma_gf.MakeRef(sigma_fes,x.GetBlock(1),0);
518 
519  int dofs = u_fes->GlobalTrueVSize()
520  + sigma_fes->GlobalTrueVSize()
521  + hatu_fes->GlobalTrueVSize()
522  + hatf_fes->GlobalTrueVSize();
523 
524  double L2Error = 0.0;
525  double rate_err = 0.0;
526  if (exact_known)
527  {
528  double u_err = u_gf.ComputeL2Error(uex);
529  double sigma_err = sigma_gf.ComputeL2Error(sigmaex);
530  L2Error = sqrt(u_err*u_err + sigma_err*sigma_err);
531  rate_err = (it) ? dim*log(err0/L2Error)/log((double)dof0/dofs) : 0.0;
532  err0 = L2Error;
533  }
534  double rate_res = (it) ? dim*log(res0/gresidual)/log((double)dof0/dofs) : 0.0;
535 
536  res0 = gresidual;
537  dof0 = dofs;
538 
539  if (myid == 0)
540  {
541  std::ios oldState(nullptr);
542  oldState.copyfmt(std::cout);
543  std::cout << std::right << std::setw(5) << it << " | "
544  << std::setw(10) << dof0 << " | ";
545  if (exact_known)
546  {
547  std::cout << std::setprecision(3) << std::setw(10)
548  << std::scientific << err0 << " | "
549  << std::setprecision(2)
550  << std::setw(6) << std::fixed << rate_err << " | " ;
551  }
552  std::cout << std::setprecision(3)
553  << std::setw(10) << std::scientific << res0 << " | "
554  << std::setprecision(2)
555  << std::setw(6) << std::fixed << rate_res << " | "
556  << std::setw(6) << std::fixed << num_iter << " | "
557  << std::endl;
558  std::cout.copyfmt(oldState);
559  }
560 
561  if (visualization)
562  {
563  const char * keys = (it == 0 && dim == 2) ? "cgRjmlk\n" : nullptr;
564  char vishost[] = "localhost";
565  int visport = 19916;
566  VisualizeField(u_out,vishost, visport, u_gf,
567  "Numerical u", 0,0, 500, 500, keys);
568  VisualizeField(sigma_out,vishost, visport, sigma_gf,
569  "Numerical flux", 501,0,500, 500, keys);
570  }
571 
572  if (paraview)
573  {
574  paraview_dc->SetCycle(it);
575  paraview_dc->SetTime((double)it);
576  paraview_dc->Save();
577  }
578 
579  if (it == ref)
580  {
581  break;
582  }
583 
584  pmesh.GeneralRefinement(elements_to_refine,1,1);
585  for (int i =0; i<trial_fes.Size(); i++)
586  {
587  trial_fes[i]->Update(false);
588  }
589  a->Update();
590 
591  coeff_fes->Update();
592  c1_gf.Update();
593  c2_gf.Update();
594  setup_test_norm_coeffs(c1_gf,c2_gf);
595  }
596 
597  if (paraview)
598  {
599  delete paraview_dc;
600  }
601 
602  delete coeff_fes;
603  delete coeff_fec;
604  delete a;
605  delete tau_fec;
606  delete v_fec;
607  delete hatf_fes;
608  delete hatf_fec;
609  delete hatu_fes;
610  delete hatu_fec;
611  delete sigma_fec;
612  delete sigma_fes;
613  delete u_fec;
614  delete u_fes;
615 
616  return 0;
617 }
618 
619 double exact_u(const Vector & X)
620 {
621  double x = X[0];
622  double y = X[1];
623  double z = 0.;
624  if (X.Size() == 3) { z = X[2]; }
625  switch (prob)
626  {
627  case sinusoidal:
628  {
629  double alpha = M_PI * (x + y + z);
630  return sin(alpha);
631  }
632  break;
633  case EJ:
634  {
635  double alpha = sqrt(1. + 4. * epsilon * epsilon * M_PI * M_PI);
636  double r1 = (1. + alpha) / (2.*epsilon);
637  double r2 = (1. - alpha) / (2.*epsilon);
638  double denom = exp(-r2) - exp(-r1);
639 
640  double g1 = exp(r2*(x-1.));
641  double g2 = exp(r1*(x-1.));
642  double g = g1-g2;
643  return g * cos(M_PI * y)/denom;
644  }
645  break;
646  case curved_streamlines:
647  {
648  double r = sqrt(x*x+y*y);
649  return atan((1.0-r)/epsilon);
650  }
651  break;
652  default:
653  MFEM_ABORT("Wrong code path");
654  return 1;
655  break;
656  }
657 }
658 
659 void exact_gradu(const Vector & X, Vector & du)
660 {
661  double x = X[0];
662  double y = X[1];
663  double z = 0.;
664  if (X.Size() == 3) { z = X[2]; }
665  du.SetSize(X.Size());
666 
667  switch (prob)
668  {
669  case sinusoidal:
670  {
671  double alpha = M_PI * (x + y + z);
672  for (int i = 0; i<du.Size(); i++)
673  {
674  du[i] = M_PI * cos(alpha);
675  }
676  }
677  break;
678  case EJ:
679  {
680  double alpha = sqrt(1. + 4. * epsilon * epsilon * M_PI * M_PI);
681  double r1 = (1. + alpha) / (2.*epsilon);
682  double r2 = (1. - alpha) / (2.*epsilon);
683  double denom = exp(-r2) - exp(-r1);
684 
685  double g1 = exp(r2*(x-1.));
686  double g1_x = r2*g1;
687  double g2 = exp(r1*(x-1.));
688  double g2_x = r1*g2;
689  double g = g1-g2;
690  double g_x = g1_x - g2_x;
691 
692  double u_x = g_x * cos(M_PI * y)/denom;
693  double u_y = -M_PI * g * sin(M_PI*y)/denom;
694  du[0] = u_x;
695  du[1] = u_y;
696  }
697  break;
698  case curved_streamlines:
699  {
700  double r = sqrt(x*x+y*y);
701  double alpha = -2.0*r + r*r + epsilon*epsilon + 1;
702  double denom = r*alpha;
703  du[0] = - x* epsilon / denom;
704  du[1] = - y* epsilon / denom;
705  }
706  break;
707  default:
708  MFEM_ABORT("Wrong code path");
709  break;
710  }
711 }
712 
713 double exact_laplacian_u(const Vector & X)
714 {
715  double x = X[0];
716  double y = X[1];
717  double z = 0.;
718  if (X.Size() == 3) { z = X[2]; }
719  switch (prob)
720  {
721  case sinusoidal:
722  {
723  double alpha = M_PI * (x + y + z);
724  double u = sin(alpha);
725  return - M_PI*M_PI * u * X.Size();
726  }
727  break;
728  case EJ:
729  {
730  double alpha = sqrt(1. + 4. * epsilon * epsilon * M_PI * M_PI);
731  double r1 = (1. + alpha) / (2.*epsilon);
732  double r2 = (1. - alpha) / (2.*epsilon);
733  double denom = exp(-r2) - exp(-r1);
734 
735  double g1 = exp(r2*(x-1.));
736  double g1_x = r2*g1;
737  double g1_xx = r2*g1_x;
738  double g2 = exp(r1*(x-1.));
739  double g2_x = r1*g2;
740  double g2_xx = r1*g2_x;
741  double g = g1-g2;
742  double g_xx = g1_xx - g2_xx;
743 
744  double u = g * cos(M_PI * y)/denom;
745  double u_xx = g_xx * cos(M_PI * y)/denom;
746  double u_yy = -M_PI * M_PI * u;
747  return u_xx + u_yy;
748  }
749  break;
750  case curved_streamlines:
751  {
752  double r = sqrt(x*x+y*y);
753  double alpha = -2.0*r + r*r + epsilon*epsilon + 1;
754  return epsilon * (r*r - epsilon*epsilon - 1.0) / (r*alpha*alpha);
755  }
756  break;
757  default:
758  MFEM_ABORT("Wrong code path");
759  return 1;
760  break;
761  }
762 }
763 
764 void exact_sigma(const Vector & X, Vector & sigma)
765 {
766  // σ = ε ∇ u
767  exact_gradu(X,sigma);
768  sigma *= epsilon;
769 }
770 
771 double exact_hatu(const Vector & X)
772 {
773  return -exact_u(X);
774 }
775 
776 void exact_hatf(const Vector & X, Vector & hatf)
777 {
778  Vector sigma;
779  Vector beta_val;
780  beta_function(X,beta_val);
781  exact_sigma(X,sigma);
782  double u = exact_u(X);
783  hatf.SetSize(X.Size());
784  for (int i = 0; i<hatf.Size(); i++)
785  {
786  hatf[i] = beta_val[i] * u - sigma[i];
787  }
788 }
789 
790 double f_exact(const Vector & X)
791 {
792  // f = - εΔu + ∇⋅(βu)
793  Vector du;
794  exact_gradu(X,du);
795  double d2u = exact_laplacian_u(X);
796 
797  Vector beta_val;
798  beta_function(X,beta_val);
799 
800  double s = 0;
801  for (int i = 0; i<du.Size(); i++)
802  {
803  s += beta_val[i] * du[i];
804  }
805  return -epsilon * d2u + s;
806 }
807 
808 double bdr_data(const Vector &X)
809 {
810  if (prob == prob_type::bdr_layer)
811  {
812  double x = X(0);
813  double y = X(1);
814 
815  if (y==0.0)
816  {
817  return -(1.0-x);
818  }
819  else if (x == 0.0)
820  {
821  return -(1.0-y);
822  }
823  else
824  {
825  return 0.0;
826  }
827  }
828  else
829  {
830  return exact_hatu(X);
831  }
832 }
833 
834 void beta_function(const Vector & X, Vector & beta_val)
835 {
836  beta_val.SetSize(2);
838  {
839  double x = X(0);
840  double y = X(1);
841  beta_val(0) = exp(x)*sin(y);
842  beta_val(1) = exp(x)*cos(y);
843  }
844  else
845  {
846  beta_val = beta;
847  }
848 }
849 
851 {
852  Array<int> vdofs;
853  ParFiniteElementSpace * fes = c1_gf.ParFESpace();
854  ParMesh * pmesh = fes->GetParMesh();
855  for (int i = 0; i < pmesh->GetNE(); i++)
856  {
857  double volume = pmesh->GetElementVolume(i);
858  double c1 = min(epsilon/volume, 1.);
859  double c2 = min(1./epsilon, 1./volume);
860  fes->GetElementDofs(i,vdofs);
861  c1_gf.SetSubVector(vdofs,c1);
862  c2_gf.SetSubVector(vdofs,c2);
863  }
864 }
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
int visport
Conjugate gradient method.
Definition: solvers.hpp:493
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1743
double exact_u(const Vector &X)
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
The Auxiliary-space Divergence Solver in hypre.
Definition: hypre.hpp:1820
void exact_sigma(const Vector &X, Vector &sigma)
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
void SetDataFormat(VTKFormat fmt)
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void beta_function(const Vector &X, Vector &beta_val)
Array< int > & RowOffsets()
Return the row offsets for block starts.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
void setup_test_norm_coeffs(ParGridFunction &c1_gf, ParGridFunction &c2_gf)
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
Helper class for ParaView visualization data.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3323
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:718
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
Class representing the parallel weak formulation. (Convenient for DPG Equations)
Definition: pweakform.hpp:25
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the ParGridFunction.
Definition: pgridfunc.cpp:96
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
STL namespace.
(Q div u, div v) for RT elements
void exact_gradu(const Vector &X, Vector &du)
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
double bdr_data(const Vector &X)
Matrix coefficient defined as the outer product of two vector coefficients.
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Definition: pgridfunc.cpp:111
bool exact_known
Definition: ex25.cpp:144
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Definition: gridfunc.cpp:2694
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
char vishost[]
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9888
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:90
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition: solvers.hpp:250
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition: fe_coll.hpp:319
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
double epsilon
void SetHighOrderOutput(bool high_order_output_)
void SetTime(double t)
Set physical time (for time-dependent simulations)
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: pgridfunc.hpp:285
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
virtual DofTransformation * GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: pfespace.cpp:470
A general vector function coefficient.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void SetRelTol(double rtol)
Definition: solvers.hpp:199
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Definition: pgridfunc.cpp:669
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Vector coefficient defined as a product of scalar and vector coefficients.
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition: array.cpp:103
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:925
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:285
double a
Definition: lissajous.cpp:41
int main(int argc, char *argv[])
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
int dim
Definition: ex24.cpp:53
double exact_laplacian_u(const Vector &X)
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:835
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:1102
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void SetLevelsOfDetail(int levels_of_detail_)
Vector beta
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:678
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition: fe_coll.hpp:434
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
double sigma(const Vector &x)
Definition: maxwell.cpp:102
virtual void Save() override
RefCoord s[3]
prob_type prob
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Class for parallel grid function.
Definition: pgridfunc.hpp:32
prob_type
Definition: ex25.cpp:146
double exact_hatu(const Vector &X)
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
A class to handle Block systems in a matrix-free implementation.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:9820
Class for parallel meshes.
Definition: pmesh.hpp:32
double GetElementVolume(int i)
Definition: mesh.cpp:119
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
void exact_hatf(const Vector &X, Vector &hatf)
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:93
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
double f(const Vector &p)
double f_exact(const Vector &X)