MFEM  v4.6.0
Finite element discretization library
pacoustics.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 acoustics example
13 //
14 // Compile with: make pacoustics
15 //
16 // sample runs
17 
18 // mpirun -np 4 pacoustics -o 3 -m ../../data/star.mesh -sref 1 -pref 2 -rnum 1.9 -sc -prob 0
19 // mpirun -np 4 pacoustics -o 3 -m ../../data/inline-quad.mesh -sref 1 -pref 2 -rnum 5.2 -sc -prob 1
20 // mpirun -np 4 pacoustics -o 4 -m ../../data/inline-tri.mesh -sref 1 -pref 2 -rnum 7.1 -sc -prob 1
21 // mpirun -np 4 pacoustics -o 2 -m ../../data/inline-hex.mesh -sref 0 -pref 1 -rnum 1.9 -sc -prob 0
22 // mpirun -np 4 pacoustics -o 3 -m ../../data/inline-quad.mesh -sref 2 -pref 1 -rnum 7.1 -sc -prob 2
23 // mpirun -np 4 pacoustics -o 2 -m ../../data/inline-hex.mesh -sref 0 -pref 1 -rnum 4.1 -sc -prob 2
24 // mpirun -np 4 pacoustics -o 3 -m meshes/scatter.mesh -sref 1 -pref 1 -rnum 7.1 -sc -prob 3
25 // mpirun -np 4 pacoustics -o 4 -m meshes/scatter.mesh -sref 1 -pref 1 -rnum 10.1 -sc -prob 4
26 // mpirun -np 4 pacoustics -o 4 -m meshes/scatter.mesh -sref 1 -pref 1 -rnum 12.1 -sc -prob 5
27 
28 // AMR runs
29 // mpirun -np 4 pacoustics -o 3 -m meshes/scatter.mesh -sref 0 -pref 10 -theta 0.75 -rnum 10.1 -sc -prob 3
30 // mpirun -np 4 pacoustics -o 3 -m meshes/scatter.mesh -sref 0 -pref 12 -theta 0.75 -rnum 20.1 -sc -prob 3
31 
32 // Description:
33 // This example code demonstrates the use of MFEM to define and solve
34 // the "ultraweak" (UW) DPG formulation for the Helmholtz problem
35 
36 // - Δ p - ω² p = f̃ , in Ω
37 // p = p₀, on ∂Ω
38 
39 // It solves the following kinds of problems
40 // 1) Known exact solutions with error convergence rates
41 // a) f̃ = 0 and p₀ is a plane wave
42 // b) A manufactured solution problem where p_exact is a Gaussian beam
43 // 2) PML problems
44 // a) Gaussian beam scattering from a square
45 // b) Plane wave scattering from a square
46 // c) Point Source
47 
48 // The DPG UW deals with the First Order System
49 // ∇ p + i ω u = 0, in Ω
50 // ∇⋅u + i ω p = f, in Ω (1)
51 // p = p₀, in ∂Ω
52 // where f:=f̃/(i ω)
53 
54 // The ultraweak-DPG formulation is obtained by integration by parts of both
55 // equations and the introduction of trace unknowns on the mesh skeleton
56 
57 // p ∈ L²(Ω), u ∈ (L²(Ω))ᵈⁱᵐ
58 // p̂ ∈ H^1/2(Ω), û ∈ H^-1/2(Ω)
59 // -(p,∇⋅v) + i ω (u,v) + <p̂,v⋅n> = 0, ∀ v ∈ H(div,Ω)
60 // -(u,∇ q) + i ω (p,q) + <û,q > = (f,q) ∀ q ∈ H^1(Ω)
61 // p̂ = p₀ on ∂Ω
62 
63 // Note:
64 // p̂ := p, û := u on the mesh skeleton
65 
66 // -------------------------------------------------------------
67 // | | p | u | p̂ | û | RHS |
68 // -------------------------------------------------------------
69 // | v | -(p, ∇⋅v) | i ω (u,v) | < p̂, v⋅n> | | |
70 // | | | | | | |
71 // | q | i ω (p,q) |-(u , ∇ q) | | < û,q > | (f,q) |
72 
73 // where (q,v) ∈ H¹(Ω) × H(div,Ω)
74 
75 // Here we use the "Adjoint Graph" norm on the test space i.e.,
76 // ||(q,v)||²ᵥ = ||A^*(q,v)||² + ||(q,v)||² where A is the
77 // acoustics operator defined by (1)
78 
79 // The PML formulation is
80 
81 // - ∇⋅(|J| J⁻¹ J⁻ᵀ ∇ p) - ω² |J| p = f
82 
83 // where J is the Jacobian of the stretching map and |J| its determinant.
84 
85 // The first order system reads
86 
87 // ∇ p + i ω α u = 0, in Ω
88 // ∇⋅u + i ω β p = f, in Ω (2)
89 // p = p₀, in ∂Ω
90 // where f:=f̃/(i ω), α:= Jᵀ J / |J|, β:= |J|
91 
92 // and the ultraweak DPG formulation
93 //
94 // p ∈ L²(Ω), u ∈ (L²(Ω))ᵈⁱᵐ
95 // p̂ ∈ H^1/2(Ω), û ∈ H^-1/2(Ω)
96 // -(p, ∇⋅v) + i ω (α u , v) + < p̂, v⋅n> = 0, ∀ v ∈ H(div,Ω)
97 // -(u , ∇ q) + i ω (β p , q) + < û, q > = (f,q) ∀ q ∈ H¹(Ω)
98 // p̂ = p₀ on ∂Ω
99 
100 // Note:
101 // p̂ := p on Γₕ (skeleton)
102 // û := u on Γₕ
103 
104 // ----------------------------------------------------------------
105 // | | p | u | p̂ | û | RHS |
106 // ----------------------------------------------------------------
107 // | v | -(p, ∇⋅v) | i ω (α u,v) | < p̂, v⋅n> | | |
108 // | | | | | | |
109 // | q | i ω (β p,q) |-(u , ∇ q) | | < û,q > | (f,q) |
110 
111 // where (q,v) ∈ H¹(Ω) × H(div,Ω)
112 
113 // Finally the test norm is defined by the adjoint operator of (2) i.e.,
114 
115 // ||(q,v)||²ᵥ = ||A^*(q,v)||² + ||(q,v)||²
116 
117 // where A is the operator defined by (2)
118 
119 // For more information see https://doi.org/10.1016/j.camwa.2017.06.044
120 
121 #include "mfem.hpp"
122 #include "util/pcomplexweakform.hpp"
123 #include "util/pml.hpp"
124 #include "../common/mfem-common.hpp"
125 #include <fstream>
126 #include <iostream>
127 
128 using namespace std;
129 using namespace mfem;
130 using namespace mfem::common;
131 
132 complex<double> acoustics_solution(const Vector & X);
133 void acoustics_solution_grad(const Vector & X,vector<complex<double>> &dp);
134 complex<double> acoustics_solution_laplacian(const Vector & X);
135 
136 double p_exact_r(const Vector &x);
137 double p_exact_i(const Vector &x);
138 void u_exact_r(const Vector &x, Vector & u);
139 void u_exact_i(const Vector &x, Vector & u);
140 double rhs_func_r(const Vector &x);
141 double rhs_func_i(const Vector &x);
142 void gradp_exact_r(const Vector &x, Vector &gradu);
143 void gradp_exact_i(const Vector &x, Vector &gradu);
144 double divu_exact_r(const Vector &x);
145 double divu_exact_i(const Vector &x);
146 double d2_exact_r(const Vector &x);
147 double d2_exact_i(const Vector &x);
148 double hatp_exact_r(const Vector & X);
149 double hatp_exact_i(const Vector & X);
150 void hatu_exact_r(const Vector & X, Vector & hatu);
151 void hatu_exact_i(const Vector & X, Vector & hatu);
152 double source_function(const Vector &x);
153 
154 int dim;
155 double omega;
156 
158 {
165 };
166 
167 static const char *enum_str[] =
168 {
169  "plane_wave",
170  "gaussian_beam",
171  "pml_general",
172  "pml_beam_scatter",
173  "pml_plane_wave_scatter",
174  "pml_pointsource"
175 };
176 
178 
179 int main(int argc, char *argv[])
180 {
181  Mpi::Init();
182  int myid = Mpi::WorldRank();
183  Hypre::Init();
184 
185  const char *mesh_file = "../../data/inline-quad.mesh";
186  int order = 1;
187  int delta_order = 1;
188  bool visualization = true;
189  double rnum=1.0;
190  double theta = 0.0;
191  bool static_cond = false;
192  int iprob = 0;
193  int sr = 0;
194  int pr = 0;
195  bool exact_known = false;
196  bool with_pml = false;
197  bool paraview = false;
198 
199  OptionsParser args(argc, argv);
200  args.AddOption(&mesh_file, "-m", "--mesh",
201  "Mesh file to use.");
202  args.AddOption(&order, "-o", "--order",
203  "Finite element order (polynomial degree)");
204  args.AddOption(&rnum, "-rnum", "--number-of-wavelengths",
205  "Number of wavelengths");
206  args.AddOption(&iprob, "-prob", "--problem", "Problem case"
207  " 0: plane wave, 1: Gaussian beam, 2: Generic PML,"
208  " 3: Scattering of a Gaussian beam"
209  " 4: Scattering of a plane wave, 5: Point source");
210  args.AddOption(&delta_order, "-do", "--delta-order",
211  "Order enrichment for DPG test space.");
212  args.AddOption(&theta, "-theta", "--theta",
213  "Theta parameter for AMR");
214  args.AddOption(&sr, "-sref", "--serial-ref",
215  "Number of parallel refinements.");
216  args.AddOption(&pr, "-pref", "--parallel-ref",
217  "Number of parallel refinements.");
218  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
219  "--no-static-condensation", "Enable static condensation.");
220  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
221  "--no-visualization",
222  "Enable or disable GLVis visualization.");
223  args.AddOption(&paraview, "-paraview", "--paraview", "-no-paraview",
224  "--no-paraview",
225  "Enable or disable ParaView visualization.");
226  args.Parse();
227  if (!args.Good())
228  {
229  if (myid == 0)
230  {
231  args.PrintUsage(cout);
232  }
233  return 1;
234  }
235 
236  if (iprob > 5) { iprob = 0; }
237  prob = (prob_type)iprob;
238  omega = 2.*M_PI*rnum;
239 
240  if (prob > 1)
241  {
242  with_pml = true;
243  if (prob > 2) { mesh_file = "meshes/scatter.mesh"; }
244  }
245  else
246  {
247  exact_known = true;
248  }
249 
250  if (myid == 0)
251  {
252  args.PrintOptions(cout);
253  }
254 
255  Mesh mesh(mesh_file, 1, 1);
256 
257  for (int i = 0; i<sr; i++)
258  {
259  mesh.UniformRefinement();
260  }
261  dim = mesh.Dimension();
262  MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example");
263 
264  CartesianPML * pml = nullptr;
265  if (with_pml)
266  {
267  Array2D<double> length(dim, 2); length = 0.125;
268  pml = new CartesianPML(&mesh,length);
269  pml->SetOmega(omega);
270  }
271 
272  mesh.EnsureNCMesh(true);
273  ParMesh pmesh(MPI_COMM_WORLD, mesh);
274  mesh.Clear();
275 
276  Array<int> attr;
277  Array<int> attrPML;
278  // PML element attribute marker
279  if (pml) { pml->SetAttributes(&pmesh, &attr, &attrPML); }
280 
281  // Define spaces
282  enum TrialSpace
283  {
284  p_space = 0,
285  u_space = 1,
286  hatp_space = 2,
287  hatu_space = 3
288  };
289  enum TestSpace
290  {
291  q_space = 0,
292  v_space = 1
293  };
294 
295  // L2 space for p
296  FiniteElementCollection *p_fec = new L2_FECollection(order-1,dim);
297  ParFiniteElementSpace *p_fes = new ParFiniteElementSpace(&pmesh,p_fec);
298 
299  // Vector L2 space for u
300  FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
301  ParFiniteElementSpace *u_fes = new ParFiniteElementSpace(&pmesh,u_fec, dim);
302 
303  // H^1/2 space for p̂
304  FiniteElementCollection * hatp_fec = new H1_Trace_FECollection(order,dim);
305  ParFiniteElementSpace *hatp_fes = new ParFiniteElementSpace(&pmesh,hatp_fec);
306 
307  // H^-1/2 space for û
308  FiniteElementCollection * hatu_fec = new RT_Trace_FECollection(order-1,dim);
309  ParFiniteElementSpace *hatu_fes = new ParFiniteElementSpace(&pmesh,hatu_fec);
310 
311  // testspace fe collections
312  int test_order = order+delta_order;
313  FiniteElementCollection * q_fec = new H1_FECollection(test_order, dim);
314  FiniteElementCollection * v_fec = new RT_FECollection(test_order-1, dim);
315 
318  trial_fes.Append(p_fes);
319  trial_fes.Append(u_fes);
320  trial_fes.Append(hatp_fes);
321  trial_fes.Append(hatu_fes);
322  test_fec.Append(q_fec);
323  test_fec.Append(v_fec);
324 
325  // Bilinear form Coefficients
326  Coefficient * omeg_cf = nullptr;
327  Coefficient * negomeg_cf = nullptr;
328  Coefficient * omeg2_cf = nullptr;
329 
330  ConstantCoefficient one(1.0);
331  ConstantCoefficient negone(-1.0);
334  ConstantCoefficient negomeg(-omega);
335 
336  if (pml)
337  {
338  omeg_cf = new RestrictedCoefficient(omeg,attr);
339  negomeg_cf = new RestrictedCoefficient(negomeg,attr);
340  omeg2_cf = new RestrictedCoefficient(omeg2,attr);
341  }
342  else
343  {
344  omeg_cf = &omeg;
345  negomeg_cf = &negomeg;
346  omeg2_cf = &omeg2;
347  }
348 
349  // PML coefficients
350  PmlCoefficient detJ_r(detJ_r_function,pml);
351  PmlCoefficient detJ_i(detJ_i_function,pml);
352  PmlCoefficient abs_detJ_2(abs_detJ_2_function,pml);
353  ProductCoefficient omeg_detJ_r(omeg,detJ_r);
354  ProductCoefficient omeg_detJ_i(omeg,detJ_i);
355  ProductCoefficient negomeg_detJ_r(negomeg,detJ_r);
356  ProductCoefficient negomeg_detJ_i(negomeg,detJ_i);
357  ProductCoefficient omeg2_abs_detJ_2(omeg2,abs_detJ_2);
358  RestrictedCoefficient omeg_detJ_r_restr(omeg_detJ_r,attrPML);
359  RestrictedCoefficient omeg_detJ_i_restr(omeg_detJ_i,attrPML);
360  RestrictedCoefficient negomeg_detJ_r_restr(negomeg_detJ_r,attrPML);
361  RestrictedCoefficient negomeg_detJ_i_restr(negomeg_detJ_i,attrPML);
362  RestrictedCoefficient omeg2_abs_detJ_2_restr(omeg2_abs_detJ_2,attrPML);
363  PmlMatrixCoefficient Jt_J_detJinv_r(dim, Jt_J_detJinv_r_function,pml);
364  PmlMatrixCoefficient Jt_J_detJinv_i(dim, Jt_J_detJinv_i_function,pml);
365  PmlMatrixCoefficient abs_Jt_J_detJinv_2(dim, abs_Jt_J_detJinv_2_function,pml);
366  ScalarMatrixProductCoefficient omeg_Jt_J_detJinv_r(omeg,Jt_J_detJinv_r);
367  ScalarMatrixProductCoefficient omeg_Jt_J_detJinv_i(omeg,Jt_J_detJinv_i);
368  ScalarMatrixProductCoefficient negomeg_Jt_J_detJinv_r(negomeg,Jt_J_detJinv_r);
369  ScalarMatrixProductCoefficient negomeg_Jt_J_detJinv_i(negomeg,Jt_J_detJinv_i);
370  ScalarMatrixProductCoefficient omeg2_abs_Jt_J_detJinv_2(omeg2,
371  abs_Jt_J_detJinv_2);
372  MatrixRestrictedCoefficient omeg_Jt_J_detJinv_r_restr(omeg_Jt_J_detJinv_r,
373  attrPML);
374  MatrixRestrictedCoefficient omeg_Jt_J_detJinv_i_restr(omeg_Jt_J_detJinv_i,
375  attrPML);
376  MatrixRestrictedCoefficient negomeg_Jt_J_detJinv_r_restr(negomeg_Jt_J_detJinv_r,
377  attrPML);
378  MatrixRestrictedCoefficient negomeg_Jt_J_detJinv_i_restr(negomeg_Jt_J_detJinv_i,
379  attrPML);
380  MatrixRestrictedCoefficient omeg2_abs_Jt_J_detJinv_2_restr(
381  omeg2_abs_Jt_J_detJinv_2,attrPML);
382 
383  ParComplexDPGWeakForm * a = new ParComplexDPGWeakForm(trial_fes,test_fec);
384  a->StoreMatrices(); // needed for AMR
385 
386  // Trial integrators
387  // Integrators not in PML
388  // i ω (p,q)
389  a->AddTrialIntegrator(nullptr,new MixedScalarMassIntegrator(*omeg_cf),
390  TrialSpace::p_space,TestSpace::q_space);
391  // -(u , ∇ q)
392  a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(negone)),
393  nullptr,TrialSpace::u_space,TestSpace::q_space);
394  // -(p, ∇⋅v)
395  a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(one),nullptr,
396  TrialSpace::p_space,TestSpace::v_space);
397  // i ω (u,v)
398  a->AddTrialIntegrator(nullptr,
399  new TransposeIntegrator(new VectorFEMassIntegrator(*omeg_cf)),
400  TrialSpace::u_space,TestSpace::v_space);
401  // < p̂, v⋅n>
402  a->AddTrialIntegrator(new NormalTraceIntegrator,nullptr,
403  TrialSpace::hatp_space,TestSpace::v_space);
404  // < û,q >
405  a->AddTrialIntegrator(new TraceIntegrator,nullptr,
406  TrialSpace::hatu_space,TestSpace::q_space);
407 
408  // test integrators
409  // (∇q,∇δq)
410  a->AddTestIntegrator(new DiffusionIntegrator(one),nullptr,
411  TestSpace::q_space, TestSpace::q_space);
412  // (q,δq)
413  a->AddTestIntegrator(new MassIntegrator(one),nullptr,
414  TestSpace::q_space, TestSpace::q_space);
415  // (∇⋅v,∇⋅δv)
416  a->AddTestIntegrator(new DivDivIntegrator(one),nullptr,
417  TestSpace::v_space, TestSpace::v_space);
418  // (v,δv)
419  a->AddTestIntegrator(new VectorFEMassIntegrator(one),nullptr,
420  TestSpace::v_space, TestSpace::v_space);
421  // -i ω (∇q,δv)
422  a->AddTestIntegrator(nullptr,new MixedVectorGradientIntegrator(*negomeg_cf),
423  TestSpace::q_space, TestSpace::v_space);
424  // i ω (v,∇ δq)
425  a->AddTestIntegrator(nullptr,
426  new MixedVectorWeakDivergenceIntegrator(*negomeg_cf),
427  TestSpace::v_space, TestSpace::q_space);
428  // ω^2 (v,δv)
429  a->AddTestIntegrator(new VectorFEMassIntegrator(*omeg2_cf),nullptr,
430  TestSpace::v_space, TestSpace::v_space);
431  // - i ω (∇⋅v,δq)
432  a->AddTestIntegrator(nullptr,new VectorFEDivergenceIntegrator(*negomeg_cf),
433  TestSpace::v_space, TestSpace::q_space);
434  // i ω (q,∇⋅v)
435  a->AddTestIntegrator(nullptr,new MixedScalarWeakGradientIntegrator(*negomeg_cf),
436  TestSpace::q_space, TestSpace::v_space);
437  // ω^2 (q,δq)
438  a->AddTestIntegrator(new MassIntegrator(*omeg2_cf),nullptr,
439  TestSpace::q_space, TestSpace::q_space);
440 
441  // integrators in the PML region
442  // Custom integration rule for the test space in the PML region
443  const IntegrationRule &ir = IntRules.Get(pmesh.GetElementGeometry(0),
444  2*test_order + 1);
445  if (pml)
446  {
447  // Trial integrators
448  // i ω (p,q) = i ω ( (β_r p,q) + i (β_i p,q) )
449  // = (- ω b_i p ) + i (ω β_r p,q)
450  a->AddTrialIntegrator(new MixedScalarMassIntegrator(negomeg_detJ_i_restr),
451  new MixedScalarMassIntegrator(omeg_detJ_r_restr),
452  TrialSpace::p_space,TestSpace::q_space);
453 
454  // i ω (α u,v) = i ω ( (α_re u,v) + i (α_im u,v) )
455  // = (-ω a_im u,v) + i (ω a_re u, v)
456  a->AddTrialIntegrator(new TransposeIntegrator(
457  new VectorFEMassIntegrator(negomeg_Jt_J_detJinv_i_restr)),
459  new VectorFEMassIntegrator(omeg_Jt_J_detJinv_r_restr)),
460  TrialSpace::u_space,TestSpace::v_space);
461  // Test integrators
462  // -i ω (α ∇q,δv) = -i ω ( (α_r ∇q,δv) + i (α_i ∇q,δv) )
463  // = (ω α_i ∇q,δv) + i (-ω α_r ∇q,δv)
465  omeg_Jt_J_detJinv_i_restr);
466  integ0_r->SetIntegrationRule(ir);
468  negomeg_Jt_J_detJinv_r_restr);
469  integ0_i->SetIntegrationRule(ir);
470  a->AddTestIntegrator(integ0_r, integ0_i,
471  TestSpace::q_space,TestSpace::v_space);
472 
473  // i ω (α^* v,∇ δq) = i ω (ᾱ v,∇ δq) (since α is diagonal)
474  // = i ω ( (α_r v,∇ δq) - i (α_i v,∇ δq)
475  // = (ω α_i v, ∇ δq) + i (ω α_r v,∇ δq )
476  a->AddTestIntegrator(new MixedVectorWeakDivergenceIntegrator(
477  negomeg_Jt_J_detJinv_i_restr),
478  new MixedVectorWeakDivergenceIntegrator(negomeg_Jt_J_detJinv_r_restr),
479  TestSpace::v_space,TestSpace::q_space);
480 
481  // ω^2 (|α|^2 v,δv) α α^* = |α|^2 since α is diagonal
483  omeg2_abs_Jt_J_detJinv_2_restr);
484  integ1->SetIntegrationRule(ir);
485  a->AddTestIntegrator(integ1, nullptr,TestSpace::v_space,TestSpace::v_space);
486 
487  // - i ω (β ∇⋅v,δq) = - i ω ( (β_re ∇⋅v,δq) + i (β_im ∇⋅v,δq) )
488  // = (ω β_im ∇⋅v,δq) + i (-ω β_re ∇⋅v,δq )
489  a->AddTestIntegrator(new VectorFEDivergenceIntegrator(omeg_detJ_i_restr),
490  new VectorFEDivergenceIntegrator(negomeg_detJ_r_restr),
491  TestSpace::v_space,TestSpace::q_space);
492 
493  // i ω (β̄ q,∇⋅v) = i ω ( (β_re ∇⋅v,δq) - i (β_im ∇⋅v,δq) )
494  // = (ω β_im ∇⋅v,δq) + i (ω β_re ∇⋅v,δq )
495  a->AddTestIntegrator(new MixedScalarWeakGradientIntegrator(
496  negomeg_detJ_i_restr),
497  new MixedScalarWeakGradientIntegrator(negomeg_detJ_r_restr),
498  TestSpace::q_space,TestSpace::v_space);
499 
500  // ω^2 (β̄ β q,δq) = (ω^2 |β|^2 )
501  MassIntegrator * integ = new MassIntegrator(omeg2_abs_detJ_2_restr);
502  integ->SetIntegrationRule(ir);
503  a->AddTestIntegrator(integ,nullptr,
504  TestSpace::q_space,TestSpace::q_space);
505  }
506 
507  // RHS
512  {
513  a->AddDomainLFIntegrator(new DomainLFIntegrator(f_rhs_r),
514  new DomainLFIntegrator(f_rhs_i),
515  TestSpace::q_space);
516  }
518  {
519  a->AddDomainLFIntegrator(new DomainLFIntegrator(f_source),nullptr,
520  TestSpace::q_space);
521  }
522 
525 
526  Array<int> elements_to_refine;
527 
528  socketstream p_out_r;
529  socketstream p_out_i;
530  if (myid == 0)
531  {
532  std::cout << "\n Ref |"
533  << " Dofs |"
534  << " ω |" ;
535  if (exact_known)
536  {
537  std::cout << " L2 Error |"
538  << " Rate |" ;
539  }
540  std::cout << " Residual |"
541  << " Rate |"
542  << " PCG it |" << endl;
543  std::cout << std::string((exact_known) ? 82 : 60,'-')
544  << endl;
545  }
546 
547  double res0 = 0.;
548  double err0 = 0.;
549  int dof0 = 0;
550 
551  ParGridFunction p_r, p_i, u_r, u_i;
552 
553  ParaViewDataCollection * paraview_dc = nullptr;
554 
555  if (paraview)
556  {
557  paraview_dc = new ParaViewDataCollection(enum_str[prob], &pmesh);
558  paraview_dc->SetPrefixPath("ParaView/Acoustics");
559  paraview_dc->SetLevelsOfDetail(order);
560  paraview_dc->SetCycle(0);
561  paraview_dc->SetDataFormat(VTKFormat::BINARY);
562  paraview_dc->SetHighOrderOutput(true);
563  paraview_dc->SetTime(0.0); // set the time
564  paraview_dc->RegisterField("p_r",&p_r);
565  paraview_dc->RegisterField("p_i",&p_i);
566  paraview_dc->RegisterField("u_r",&u_r);
567  paraview_dc->RegisterField("u_i",&u_i);
568  }
569 
570  if (static_cond) { a->EnableStaticCondensation(); }
571  for (int it = 0; it<=pr; it++)
572  {
573  a->Assemble();
574 
575  Array<int> ess_tdof_list;
576  Array<int> ess_bdr;
577  if (pmesh.bdr_attributes.Size())
578  {
579  ess_bdr.SetSize(pmesh.bdr_attributes.Max());
580  ess_bdr = 1;
581  hatp_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
582  if (pml && prob>2)
583  {
584  ess_bdr = 0;
585  ess_bdr[1] = 1;
586  }
587  }
588 
589  // shift the ess_tdofs
590  for (int j = 0; j < ess_tdof_list.Size(); j++)
591  {
592  ess_tdof_list[j] += p_fes->GetTrueVSize() + u_fes->GetTrueVSize();
593  }
594 
595  Array<int> offsets(5);
596  offsets[0] = 0;
597  offsets[1] = p_fes->GetVSize();
598  offsets[2] = u_fes->GetVSize();
599  offsets[3] = hatp_fes->GetVSize();
600  offsets[4] = hatu_fes->GetVSize();
601  offsets.PartialSum();
602 
603  Vector x(2*offsets.Last());
604  x = 0.;
605 
606  if (prob!=2)
607  {
608  ParGridFunction hatp_gf_r(hatp_fes, x, offsets[2]);
609  ParGridFunction hatp_gf_i(hatp_fes, x, offsets.Last()+ offsets[2]);
610  hatp_gf_r.ProjectBdrCoefficient(hatpex_r, ess_bdr);
611  hatp_gf_i.ProjectBdrCoefficient(hatpex_i, ess_bdr);
612  }
613 
614  OperatorPtr Ah;
615  Vector X,B;
616  a->FormLinearSystem(ess_tdof_list,x,Ah, X,B);
617 
618  ComplexOperator * Ahc = Ah.As<ComplexOperator>();
619  BlockOperator * BlockA_r = dynamic_cast<BlockOperator *>(&Ahc->real());
620  BlockOperator * BlockA_i = dynamic_cast<BlockOperator *>(&Ahc->imag());
621 
622  int num_blocks = BlockA_r->NumRowBlocks();
623  Array<int> tdof_offsets(2*num_blocks+1);
624 
625  tdof_offsets[0] = 0;
626  int skip = (static_cond) ? 0 : 2;
627  int k = (static_cond) ? 2 : 0;
628  for (int i=0; i<num_blocks; i++)
629  {
630  tdof_offsets[i+1] = trial_fes[i+k]->GetTrueVSize();
631  tdof_offsets[num_blocks+i+1] = trial_fes[i+k]->GetTrueVSize();
632  }
633  tdof_offsets.PartialSum();
634 
635  BlockOperator blockA(tdof_offsets);
636  for (int i = 0; i<num_blocks; i++)
637  {
638  for (int j = 0; j<num_blocks; j++)
639  {
640  blockA.SetBlock(i,j,&BlockA_r->GetBlock(i,j));
641  blockA.SetBlock(i,j+num_blocks,&BlockA_i->GetBlock(i,j), -1.0);
642  blockA.SetBlock(i+num_blocks,j+num_blocks,&BlockA_r->GetBlock(i,j));
643  blockA.SetBlock(i+num_blocks,j,&BlockA_i->GetBlock(i,j));
644  }
645  }
646 
647  X = 0.;
648  BlockDiagonalPreconditioner M(tdof_offsets);
649  M.owns_blocks=0;
650 
651  if (!static_cond)
652  {
653  HypreBoomerAMG * solver_p = new HypreBoomerAMG((HypreParMatrix &)
654  BlockA_r->GetBlock(0,0));
655  solver_p->SetPrintLevel(0);
656  solver_p->SetSystemsOptions(dim);
657  HypreBoomerAMG * solver_u = new HypreBoomerAMG((HypreParMatrix &)
658  BlockA_r->GetBlock(1,1));
659  solver_u->SetPrintLevel(0);
660  solver_u->SetSystemsOptions(dim);
661  M.SetDiagonalBlock(0,solver_p);
662  M.SetDiagonalBlock(1,solver_u);
663  M.SetDiagonalBlock(num_blocks,solver_p);
664  M.SetDiagonalBlock(num_blocks+1,solver_u);
665  }
666 
667  HypreBoomerAMG * solver_hatp = new HypreBoomerAMG((HypreParMatrix &)
668  BlockA_r->GetBlock(skip,skip));
669  solver_hatp->SetPrintLevel(0);
670 
671  HypreSolver * solver_hatu = nullptr;
672  if (dim == 2)
673  {
674  // AMS preconditioner for 2D H(div) (trace) space
675  solver_hatu = new HypreAMS((HypreParMatrix &)BlockA_r->GetBlock(skip+1,skip+1),
676  hatu_fes);
677  dynamic_cast<HypreAMS*>(solver_hatu)->SetPrintLevel(0);
678  }
679  else
680  {
681  // ADS preconditioner for 3D H(div) (trace) space
682  solver_hatu = new HypreADS((HypreParMatrix &)BlockA_r->GetBlock(skip+1,skip+1),
683  hatu_fes);
684  dynamic_cast<HypreADS*>(solver_hatu)->SetPrintLevel(0);
685  }
686 
687  M.SetDiagonalBlock(skip,solver_hatp);
688  M.SetDiagonalBlock(skip+1,solver_hatu);
689  M.SetDiagonalBlock(skip+num_blocks,solver_hatp);
690  M.SetDiagonalBlock(skip+num_blocks+1,solver_hatu);
691 
692  CGSolver cg(MPI_COMM_WORLD);
693  cg.SetRelTol(1e-6);
694  cg.SetMaxIter(10000);
695  cg.SetPrintLevel(0);
696  cg.SetPreconditioner(M);
697  cg.SetOperator(blockA);
698  cg.Mult(B, X);
699 
700  for (int i = 0; i<num_blocks; i++)
701  {
702  delete &M.GetDiagonalBlock(i);
703  }
704 
705  int num_iter = cg.GetNumIterations();
706 
707  a->RecoverFEMSolution(X,x);
708 
709  Vector & residuals = a->ComputeResidual(x);
710 
711  double residual = residuals.Norml2();
712  double maxresidual = residuals.Max();
713  double globalresidual = residual * residual;
714  MPI_Allreduce(MPI_IN_PLACE,&maxresidual,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
715  MPI_Allreduce(MPI_IN_PLACE,&globalresidual,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
716 
717  globalresidual = sqrt(globalresidual);
718 
719  p_r.MakeRef(p_fes, x, 0);
720  p_i.MakeRef(p_fes, x, offsets.Last());
721 
722  u_r.MakeRef(u_fes,x, offsets[1]);
723  u_i.MakeRef(u_fes,x, offsets.Last()+offsets[1]);
724 
725  int dofs = 0;
726  for (int i = 0; i<trial_fes.Size(); i++)
727  {
728  dofs += trial_fes[i]->GlobalTrueVSize();
729  }
730 
731  double L2Error = 0.0;
732  double rate_err = 0.0;
733  if (exact_known)
734  {
737  double p_err_r = p_r.ComputeL2Error(p_ex_r);
738  double p_err_i = p_i.ComputeL2Error(p_ex_i);
739 
740  // Error in velocity
743 
744  double u_err_r = u_r.ComputeL2Error(u_ex_r);
745  double u_err_i = u_i.ComputeL2Error(u_ex_i);
746 
747  L2Error = sqrt(p_err_r*p_err_r + p_err_i*p_err_i
748  +u_err_r*u_err_r + u_err_i*u_err_i);
749 
750  rate_err = (it) ? dim*log(err0/L2Error)/log((double)dof0/dofs) : 0.0;
751  err0 = L2Error;
752  }
753 
754  double rate_res = (it) ? dim*log(res0/globalresidual)/log((
755  double)dof0/dofs) : 0.0;
756 
757  res0 = globalresidual;
758  dof0 = dofs;
759 
760  if (myid == 0)
761  {
762  std::ios oldState(nullptr);
763  oldState.copyfmt(std::cout);
764  std::cout << std::right << std::setw(5) << it << " | "
765  << std::setw(10) << dof0 << " | "
766  << std::setprecision(1) << std::fixed
767  << std::setw(4) << 2*rnum << " π | ";
768  if (exact_known)
769  {
770  std::cout << std::setprecision(3) << std::setw(10)
771  << std::scientific << err0 << " | "
772  << std::setprecision(2)
773  << std::setw(6) << std::fixed << rate_err << " | " ;
774  }
775  std::cout << std::setprecision(3)
776  << std::setw(10) << std::scientific << res0 << " | "
777  << std::setprecision(2)
778  << std::setw(6) << std::fixed << rate_res << " | "
779  << std::setw(6) << std::fixed << num_iter << " | "
780  << std::endl;
781  std::cout.copyfmt(oldState);
782  }
783 
784  if (visualization)
785  {
786  const char * keys = (it == 0 && dim == 2) ? "jRcml\n" : nullptr;
787  char vishost[] = "localhost";
788  int visport = 19916;
789  VisualizeField(p_out_r,vishost, visport, p_r,
790  "Numerical presure (real part)", 0, 0, 500, 500, keys);
791  VisualizeField(p_out_i,vishost, visport, p_i,
792  "Numerical presure (imaginary part)", 501, 0, 500, 500, keys);
793  }
794 
795  if (paraview)
796  {
797  paraview_dc->SetCycle(it);
798  paraview_dc->SetTime((double)it);
799  paraview_dc->Save();
800  }
801 
802  if (it == pr)
803  {
804  break;
805  }
806 
807  if (theta > 0.0)
808  {
809  elements_to_refine.SetSize(0);
810  for (int iel = 0; iel<pmesh.GetNE(); iel++)
811  {
812  if (residuals[iel] > theta * maxresidual)
813  {
814  elements_to_refine.Append(iel);
815  }
816  }
817  pmesh.GeneralRefinement(elements_to_refine,1,1);
818  }
819  else
820  {
821  pmesh.UniformRefinement();
822  }
823  if (pml) { pml->SetAttributes(&pmesh); }
824  for (int i =0; i<trial_fes.Size(); i++)
825  {
826  trial_fes[i]->Update(false);
827  }
828  a->Update();
829  }
830 
831  if (paraview)
832  {
833  delete paraview_dc;
834  }
835 
836  if (pml)
837  {
838  delete omeg_cf;
839  delete omeg2_cf;
840  delete negomeg_cf;
841  delete pml;
842  }
843  delete a;
844  delete q_fec;
845  delete v_fec;
846  delete hatp_fes;
847  delete hatp_fec;
848  delete hatu_fes;
849  delete hatu_fec;
850  delete u_fec;
851  delete p_fec;
852  delete u_fes;
853  delete p_fes;
854 
855  return 0;
856 }
857 
858 double p_exact_r(const Vector &x)
859 {
860  return acoustics_solution(x).real();
861 }
862 
863 double p_exact_i(const Vector &x)
864 {
865  return acoustics_solution(x).imag();
866 }
867 
868 double hatp_exact_r(const Vector & X)
869 {
870  return p_exact_r(X);
871 }
872 
873 double hatp_exact_i(const Vector & X)
874 {
875  return p_exact_i(X);
876 }
877 
878 void gradp_exact_r(const Vector &x, Vector &grad_r)
879 {
880  grad_r.SetSize(x.Size());
881  vector<complex<double>> grad;
882  acoustics_solution_grad(x,grad);
883  for (unsigned i = 0; i < grad.size(); i++)
884  {
885  grad_r[i] = grad[i].real();
886  }
887 }
888 
889 void gradp_exact_i(const Vector &x, Vector &grad_i)
890 {
891  grad_i.SetSize(x.Size());
892  vector<complex<double>> grad;
893  acoustics_solution_grad(x,grad);
894  for (unsigned i = 0; i < grad.size(); i++)
895  {
896  grad_i[i] = grad[i].imag();
897  }
898 }
899 
900 double d2_exact_r(const Vector &x)
901 {
902  return acoustics_solution_laplacian(x).real();
903 }
904 
905 double d2_exact_i(const Vector &x)
906 {
907  return acoustics_solution_laplacian(x).imag();
908 }
909 
910 // u = - ∇ p / (i ω )
911 // = i (∇ p_r + i * ∇ p_i) / ω
912 // = - ∇ p_i / ω + i ∇ p_r / ω
913 void u_exact_r(const Vector &x, Vector & u)
914 {
915  gradp_exact_i(x,u);
916  u *= -1./omega;
917 }
918 
919 void u_exact_i(const Vector &x, Vector & u)
920 {
921  gradp_exact_r(x,u);
922  u *= 1./omega;
923 }
924 
925 void hatu_exact_r(const Vector & X, Vector & hatu)
926 {
927  u_exact_r(X,hatu);
928 }
929 void hatu_exact_i(const Vector & X, Vector & hatu)
930 {
931  u_exact_i(X,hatu);
932 }
933 
934 // ∇⋅u = i Δ p / ω
935 // = i (Δ p_r + i * Δ p_i) / ω
936 // = - Δ p_i / ω + i Δ p_r / ω
937 
938 double divu_exact_r(const Vector &x)
939 {
940  return -d2_exact_i(x)/omega;
941 }
942 
943 double divu_exact_i(const Vector &x)
944 {
945  return d2_exact_r(x)/omega;
946 }
947 
948 // f = ∇⋅u + i ω p
949 // f_r = ∇⋅u_r - ω p_i
950 double rhs_func_r(const Vector &x)
951 {
952  double p = p_exact_i(x);
953  double divu = divu_exact_r(x);
954  return divu - omega * p;
955 }
956 
957 // f_i = ∇⋅u_i + ω p_r
958 double rhs_func_i(const Vector &x)
959 {
960  double p = p_exact_r(x);
961  double divu = divu_exact_i(x);
962  return divu + omega * p;
963 }
964 
965 complex<double> acoustics_solution(const Vector & X)
966 {
967  complex<double> zi = complex<double>(0., 1.);
968  switch (prob)
969  {
971  case plane_wave:
972  {
973  double beta = omega/std::sqrt((double)X.Size());
974  complex<double> alpha = beta * zi * X.Sum();
975  return exp(alpha);
976  }
977  break;
978  case gaussian_beam:
979  case pml_beam_scatter:
980  {
981  double rk = omega;
982  double degrees = 45;
983  double alpha = (180+degrees) * M_PI/180.;
984  double sina = sin(alpha);
985  double cosa = cos(alpha);
986  // shift the origin
987  double shift = 0.1;
988  double xprim=X(0) + shift;
989  double yprim=X(1) + shift;
990 
991  double x = xprim*sina - yprim*cosa;
992  double y = xprim*cosa + yprim*sina;
993  //wavelength
994  double rl = 2.*M_PI/rk;
995 
996  // beam waist radius
997  double w0 = 0.05;
998 
999  // function w
1000  double fact = rl/M_PI/(w0*w0);
1001  double aux = 1. + (fact*y)*(fact*y);
1002 
1003  double w = w0*sqrt(aux);
1004 
1005  double phi0 = atan(fact*y);
1006 
1007  double r = y + 1./y/(fact*fact);
1008 
1009  // pressure
1010  complex<double> ze = - x*x/(w*w) - zi*rk*y - zi * M_PI * x * x/rl/r +
1011  zi*phi0/2.;
1012  double pf = pow(2.0/M_PI/(w*w),0.25);
1013 
1014  return pf*exp(ze);
1015  }
1016  break;
1017  case pml_pointsource:
1018  {
1019  double x = X(0)-0.5;
1020  double y = X(1)-0.5;
1021  double r = sqrt(x*x + y*y);
1022  double beta = omega * r;
1023  complex<double> Ho = jn(0, beta) + zi * yn(0, beta);
1024  return 0.25*zi*Ho;
1025  }
1026  break;
1027  default:
1028  MFEM_ABORT("Should be unreachable");
1029  return 1;
1030  break;
1031  }
1032 }
1033 
1034 void acoustics_solution_grad(const Vector & X, vector<complex<double>> & dp)
1035 {
1036  dp.resize(X.Size());
1037  complex<double> zi = complex<double>(0., 1.);
1038  // initialize
1039  for (int i = 0; i<X.Size(); i++) { dp[i] = 0.0; }
1040  switch (prob)
1041  {
1043  case plane_wave:
1044  {
1045  double beta = omega/std::sqrt((double)X.Size());
1046  complex<double> alpha = beta * zi * X.Sum();
1047  complex<double> p = exp(alpha);
1048  for (int i = 0; i<X.Size(); i++)
1049  {
1050  dp[i] = zi * beta * p;
1051  }
1052  }
1053  break;
1054  case gaussian_beam:
1055  case pml_beam_scatter:
1056  {
1057  double rk = omega;
1058  double degrees = 45;
1059  double alpha = (180+degrees) * M_PI/180.;
1060  double sina = sin(alpha);
1061  double cosa = cos(alpha);
1062  // shift the origin
1063  double shift = 0.1;
1064  double xprim=X(0) + shift;
1065  double yprim=X(1) + shift;
1066 
1067  double x = xprim*sina - yprim*cosa;
1068  double y = xprim*cosa + yprim*sina;
1069  double dxdxprim = sina, dxdyprim = -cosa;
1070  double dydxprim = cosa, dydyprim = sina;
1071  //wavelength
1072  double rl = 2.*M_PI/rk;
1073 
1074  // beam waist radius
1075  double w0 = 0.05;
1076 
1077  // function w
1078  double fact = rl/M_PI/(w0*w0);
1079  double aux = 1. + (fact*y)*(fact*y);
1080 
1081  double w = w0*sqrt(aux);
1082  double dwdy = w0*fact*fact*y/sqrt(aux);
1083 
1084  double phi0 = atan(fact*y);
1085  double dphi0dy = cos(phi0)*cos(phi0)*fact;
1086 
1087  double r = y + 1./y/(fact*fact);
1088  double drdy = 1. - 1./(y*y)/(fact*fact);
1089 
1090  // pressure
1091  complex<double> ze = - x*x/(w*w) - zi*rk*y - zi * M_PI * x * x/rl/r +
1092  zi*phi0/2.;
1093 
1094  complex<double> zdedx = -2.*x/(w*w) - 2.*zi*M_PI*x/rl/r;
1095  complex<double> zdedy = 2.*x*x/(w*w*w)*dwdy - zi*rk + zi*M_PI*x*x/rl/
1096  (r*r)*drdy + zi*dphi0dy/2.;
1097 
1098  double pf = pow(2.0/M_PI/(w*w),0.25);
1099  double dpfdy = -pow(2./M_PI/(w*w),-0.75)/M_PI/(w*w*w)*dwdy;
1100 
1101  complex<double> zp = pf*exp(ze);
1102  complex<double> zdpdx = zp*zdedx;
1103  complex<double> zdpdy = dpfdy*exp(ze)+zp*zdedy;
1104 
1105  dp[0] = (zdpdx*dxdxprim + zdpdy*dydxprim);
1106  dp[1] = (zdpdx*dxdyprim + zdpdy*dydyprim);
1107  }
1108  break;
1109  default:
1110  MFEM_ABORT("Should be unreachable");
1111  break;
1112  }
1113 }
1114 
1115 complex<double> acoustics_solution_laplacian(const Vector & X)
1116 {
1117  complex<double> zi = complex<double>(0., 1.);
1118  switch (prob)
1119  {
1121  case plane_wave:
1122  {
1123  double beta = omega/std::sqrt((double)X.Size());
1124  complex<double> alpha = beta * zi * X.Sum();
1125  return dim * beta * beta * exp(alpha);
1126  }
1127  break;
1128  case gaussian_beam:
1129  case pml_beam_scatter:
1130  {
1131  double rk = omega;
1132  double degrees = 45;
1133  double alpha = (180+degrees) * M_PI/180.;
1134  double sina = sin(alpha);
1135  double cosa = cos(alpha);
1136  // shift the origin
1137  double shift = 0.1;
1138  double xprim=X(0) + shift;
1139  double yprim=X(1) + shift;
1140 
1141  double x = xprim*sina - yprim*cosa;
1142  double y = xprim*cosa + yprim*sina;
1143  double dxdxprim = sina, dxdyprim = -cosa;
1144  double dydxprim = cosa, dydyprim = sina;
1145  //wavelength
1146  double rl = 2.*M_PI/rk;
1147 
1148  // beam waist radius
1149  double w0 = 0.05;
1150 
1151  // function w
1152  double fact = rl/M_PI/(w0*w0);
1153  double aux = 1. + (fact*y)*(fact*y);
1154 
1155  double w = w0*sqrt(aux);
1156  double dwdy = w0*fact*fact*y/sqrt(aux);
1157  double d2wdydy = w0*fact*fact*(1. - (fact*y)*(fact*y)/aux)/sqrt(aux);
1158 
1159  double phi0 = atan(fact*y);
1160  double dphi0dy = cos(phi0)*cos(phi0)*fact;
1161  double d2phi0dydy = -2.*cos(phi0)*sin(phi0)*fact*dphi0dy;
1162 
1163  double r = y + 1./y/(fact*fact);
1164  double drdy = 1. - 1./(y*y)/(fact*fact);
1165  double d2rdydy = 2./(y*y*y)/(fact*fact);
1166 
1167  // pressure
1168  complex<double> ze = - x*x/(w*w) - zi*rk*y - zi * M_PI * x * x/rl/r +
1169  zi*phi0/2.;
1170 
1171  complex<double> zdedx = -2.*x/(w*w) - 2.*zi*M_PI*x/rl/r;
1172  complex<double> zdedy = 2.*x*x/(w*w*w)*dwdy - zi*rk + zi*M_PI*x*x/rl/
1173  (r*r)*drdy + zi*dphi0dy/2.;
1174  complex<double> zd2edxdx = -2./(w*w) - 2.*zi*M_PI/rl/r;
1175  complex<double> zd2edxdy = 4.*x/(w*w*w)*dwdy + 2.*zi*M_PI*x/rl/(r*r)*drdy;
1176  complex<double> zd2edydx = zd2edxdy;
1177  complex<double> zd2edydy = -6.*x*x/(w*w*w*w)*dwdy*dwdy + 2.*x*x/
1178  (w*w*w)*d2wdydy - 2.*zi*M_PI*x*x/rl/(r*r*r)*drdy*drdy
1179  + zi*M_PI*x*x/rl/(r*r)*d2rdydy + zi/2.*d2phi0dydy;
1180 
1181  double pf = pow(2.0/M_PI/(w*w),0.25);
1182  double dpfdy = -pow(2./M_PI/(w*w),-0.75)/M_PI/(w*w*w)*dwdy;
1183  double d2pfdydy = -1./M_PI*pow(2./M_PI,-0.75)*(-1.5*pow(w,-2.5)
1184  *dwdy*dwdy + pow(w,-1.5)*d2wdydy);
1185 
1186 
1187  complex<double> zp = pf*exp(ze);
1188  complex<double> zdpdx = zp*zdedx;
1189  complex<double> zdpdy = dpfdy*exp(ze)+zp*zdedy;
1190  complex<double> zd2pdxdx = zdpdx*zdedx + zp*zd2edxdx;
1191  complex<double> zd2pdxdy = zdpdy*zdedx + zp*zd2edxdy;
1192  complex<double> zd2pdydx = dpfdy*exp(ze)*zdedx + zdpdx*zdedy + zp*zd2edydx;
1193  complex<double> zd2pdydy = d2pfdydy*exp(ze) + dpfdy*exp(
1194  ze)*zdedy + zdpdy*zdedy + zp*zd2edydy;
1195 
1196 
1197  return (zd2pdxdx*dxdxprim + zd2pdydx*dydxprim)*dxdxprim
1198  + (zd2pdxdy*dxdxprim + zd2pdydy*dydxprim)*dydxprim
1199  + (zd2pdxdx*dxdyprim + zd2pdydx*dydyprim)*dxdyprim
1200  + (zd2pdxdy*dxdyprim + zd2pdydy*dydyprim)*dydyprim;
1201  }
1202  break;
1203  default:
1204  MFEM_ABORT("Should be unreachable");
1205  return 1;
1206  break;
1207  }
1208 }
1209 
1210 double source_function(const Vector &x)
1211 {
1212  Vector center(dim);
1213  center = 0.5;
1214  double r = 0.0;
1215  for (int i = 0; i < dim; ++i)
1216  {
1217  r += pow(x[i] - center[i], 2.);
1218  }
1219  double n = 5.0 * omega / M_PI;
1220  double coeff = pow(n, 2) / M_PI;
1221  double alpha = -pow(n, 2) * r;
1222  return -omega * coeff * exp(alpha)/omega;
1223 }
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
virtual Operator & real()
Real or imaginary part accessor methods.
double omega
Definition: pacoustics.cpp:155
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1743
double source_function(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
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
void SetDataFormat(VTKFormat fmt)
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
double rhs_func_r(const Vector &x)
Definition: pacoustics.cpp:950
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
double d2_exact_r(const Vector &x)
Definition: pacoustics.cpp:900
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
void SetIntegrationRule(const IntegrationRule &ir)
Prescribe a fixed IntegrationRule to use.
Definition: nonlininteg.hpp:68
Mimic the action of a complex operator using two real operators.
Helper class for ParaView visualization data.
int NumRowBlocks() const
Return the number of row blocks.
void hatu_exact_i(const Vector &X, Vector &hatu)
Definition: pacoustics.cpp:929
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
double divu_exact_r(const Vector &x)
Definition: pacoustics.cpp:938
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
complex< double > acoustics_solution_laplacian(const Vector &X)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void gradp_exact_i(const Vector &x, Vector &gradu)
Definition: pacoustics.cpp:889
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
Vector beta
void Jt_J_detJinv_i_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition: pml.cpp:203
Matrix coefficient defined as a product of a scalar coefficient and a matrix coefficient.
STL namespace.
void SetOmega(double omega_)
Definition: pml.hpp:64
complex< double > acoustics_solution(const Vector &X)
Definition: pacoustics.cpp:965
(Q div u, div v) for RT elements
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
int dim
Definition: pacoustics.cpp:154
double detJ_i_function(const Vector &x, CartesianPML *pml)
Definition: pml.cpp:166
Operator & GetDiagonalBlock(int iblock)
Return a reference to block i,i.
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
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
Derived coefficient that takes the value of the parent coefficient for the active attributes and is z...
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition: solvers.hpp:250
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
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
prob_type
Definition: pacoustics.cpp:157
void abs_Jt_J_detJinv_2_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition: pml.cpp:219
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:1286
void SetAttributes(Mesh *mesh_, Array< int > *attrNonPML=nullptr, Array< int > *attrPML=nullptr)
Mark element in the PML region.
Definition: pml.cpp:68
double abs_detJ_2_function(const Vector &x, CartesianPML *pml)
Definition: pml.cpp:176
void u_exact_r(const Vector &x, Vector &u)
Definition: pacoustics.cpp:913
Class representing the parallel weak formulation. (Convenient for DPG Equations)
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
void gradp_exact_r(const Vector &x, Vector &gradu)
Definition: pacoustics.cpp:878
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
double d2_exact_i(const Vector &x)
Definition: pacoustics.cpp:905
A general vector function coefficient.
void u_exact_i(const Vector &x, Vector &u)
Definition: pacoustics.cpp:919
Geometry::Type GetElementGeometry(int i) const
Definition: mesh.hpp:1228
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
double hatp_exact_r(const Vector &X)
Definition: pacoustics.cpp:868
double hatp_exact_i(const Vector &X)
Definition: pacoustics.cpp:873
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Definition: pgridfunc.cpp:669
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
double divu_exact_i(const Vector &x)
Definition: pacoustics.cpp:943
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
virtual Operator & imag()
double p_exact_i(const Vector &x)
Definition: pacoustics.cpp:863
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:925
Derived matrix coefficient that has the value of the parent matrix coefficient where it is active and...
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
double detJ_r_function(const Vector &x, CartesianPML *pml)
PML stretching functions: See https://doi.org/10.1006/jcph.1994.1159.
Definition: pml.cpp:156
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:285
double a
Definition: lissajous.cpp:41
Class for setting up a simple Cartesian PML region.
Definition: pml.hpp:18
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
double rhs_func_i(const Vector &x)
Definition: pacoustics.cpp:958
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
void Jt_J_detJinv_r_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition: pml.cpp:187
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
T & Last()
Return the last element in the array.
Definition: array.hpp:792
void SetLevelsOfDetail(int levels_of_detail_)
void acoustics_solution_grad(const Vector &X, vector< complex< double >> &dp)
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
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition: hypre.cpp:4957
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
void hatu_exact_r(const Vector &X, Vector &hatu)
Definition: pacoustics.cpp:925
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
virtual void Save() override
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
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 p_exact_r(const Vector &x)
Definition: pacoustics.cpp:858
prob_type prob
Definition: pacoustics.cpp:177
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
int main(int argc, char *argv[])
Definition: pacoustics.cpp:179
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327