MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pacoustics.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 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 7 -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"
123#include "util/pml.hpp"
125#include <fstream>
126#include <iostream>
127
128using namespace std;
129using namespace mfem;
130using namespace mfem::common;
131
132complex<real_t> acoustics_solution(const Vector & X);
133void acoustics_solution_grad(const Vector & X,vector<complex<real_t>> &dp);
134complex<real_t> acoustics_solution_laplacian(const Vector & X);
135
136real_t p_exact_r(const Vector &x);
137real_t p_exact_i(const Vector &x);
138void u_exact_r(const Vector &x, Vector & u);
139void u_exact_i(const Vector &x, Vector & u);
140real_t rhs_func_r(const Vector &x);
141real_t rhs_func_i(const Vector &x);
142void gradp_exact_r(const Vector &x, Vector &gradu);
143void gradp_exact_i(const Vector &x, Vector &gradu);
144real_t divu_exact_r(const Vector &x);
145real_t divu_exact_i(const Vector &x);
146real_t d2_exact_r(const Vector &x);
147real_t d2_exact_i(const Vector &x);
148real_t hatp_exact_r(const Vector & X);
149real_t hatp_exact_i(const Vector & X);
150void hatu_exact_r(const Vector & X, Vector & hatu);
151void hatu_exact_i(const Vector & X, Vector & hatu);
153
154int dim;
156
166
167static 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
179int 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 real_t rnum=1.0;
190 real_t 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<real_t> 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
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);
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,
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,
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
511 if (prob == prob_type::gaussian_beam)
512 {
513 a->AddDomainLFIntegrator(new DomainLFIntegrator(f_rhs_r),
514 new DomainLFIntegrator(f_rhs_i),
515 TestSpace::q_space);
516 }
517 if (prob == prob_type::pml_general)
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 real_t res0 = 0.;
548 real_t 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 real_t residual = residuals.Norml2();
712 real_t maxresidual = residuals.Max();
713 real_t globalresidual = residual * residual;
714 MPI_Allreduce(MPI_IN_PLACE, &maxresidual, 1, MPITypeMap<real_t>::mpi_type,
715 MPI_MAX,MPI_COMM_WORLD);
716 MPI_Allreduce(MPI_IN_PLACE, &globalresidual, 1,
717 MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
718
719 globalresidual = sqrt(globalresidual);
720
721 p_r.MakeRef(p_fes, x, 0);
722 p_i.MakeRef(p_fes, x, offsets.Last());
723
724 u_r.MakeRef(u_fes,x, offsets[1]);
725 u_i.MakeRef(u_fes,x, offsets.Last()+offsets[1]);
726
727 int dofs = 0;
728 for (int i = 0; i<trial_fes.Size(); i++)
729 {
730 dofs += trial_fes[i]->GlobalTrueVSize();
731 }
732
733 real_t L2Error = 0.0;
734 real_t rate_err = 0.0;
735 if (exact_known)
736 {
739 real_t p_err_r = p_r.ComputeL2Error(p_ex_r);
740 real_t p_err_i = p_i.ComputeL2Error(p_ex_i);
741
742 // Error in velocity
745
746 real_t u_err_r = u_r.ComputeL2Error(u_ex_r);
747 real_t u_err_i = u_i.ComputeL2Error(u_ex_i);
748
749 L2Error = sqrt(p_err_r*p_err_r + p_err_i*p_err_i
750 +u_err_r*u_err_r + u_err_i*u_err_i);
751
752 rate_err = (it) ? dim*log(err0/L2Error)/log((real_t)dof0/dofs) : 0.0;
753 err0 = L2Error;
754 }
755
756 real_t rate_res = (it) ? dim*log(res0/globalresidual)/log((
757 real_t)dof0/dofs) : 0.0;
758
759 res0 = globalresidual;
760 dof0 = dofs;
761
762 if (myid == 0)
763 {
764 std::ios oldState(nullptr);
765 oldState.copyfmt(std::cout);
766 std::cout << std::right << std::setw(5) << it << " | "
767 << std::setw(10) << dof0 << " | "
768 << std::setprecision(1) << std::fixed
769 << std::setw(4) << 2*rnum << " π | ";
770 if (exact_known)
771 {
772 std::cout << std::setprecision(3) << std::setw(10)
773 << std::scientific << err0 << " | "
774 << std::setprecision(2)
775 << std::setw(6) << std::fixed << rate_err << " | " ;
776 }
777 std::cout << std::setprecision(3)
778 << std::setw(10) << std::scientific << res0 << " | "
779 << std::setprecision(2)
780 << std::setw(6) << std::fixed << rate_res << " | "
781 << std::setw(6) << std::fixed << num_iter << " | "
782 << std::endl;
783 std::cout.copyfmt(oldState);
784 }
785
786 if (visualization)
787 {
788 const char * keys = (it == 0 && dim == 2) ? "jRcml\n" : nullptr;
789 char vishost[] = "localhost";
790 int visport = 19916;
791 VisualizeField(p_out_r,vishost, visport, p_r,
792 "Numerical presure (real part)", 0, 0, 500, 500, keys);
793 VisualizeField(p_out_i,vishost, visport, p_i,
794 "Numerical presure (imaginary part)", 501, 0, 500, 500, keys);
795 }
796
797 if (paraview)
798 {
799 paraview_dc->SetCycle(it);
800 paraview_dc->SetTime((real_t)it);
801 paraview_dc->Save();
802 }
803
804 if (it == pr)
805 {
806 break;
807 }
808
809 if (theta > 0.0)
810 {
811 elements_to_refine.SetSize(0);
812 for (int iel = 0; iel<pmesh.GetNE(); iel++)
813 {
814 if (residuals[iel] > theta * maxresidual)
815 {
816 elements_to_refine.Append(iel);
817 }
818 }
819 pmesh.GeneralRefinement(elements_to_refine,1,1);
820 }
821 else
822 {
823 pmesh.UniformRefinement();
824 }
825 if (pml) { pml->SetAttributes(&pmesh); }
826 for (int i =0; i<trial_fes.Size(); i++)
827 {
828 trial_fes[i]->Update(false);
829 }
830 a->Update();
831 }
832
833 if (paraview)
834 {
835 delete paraview_dc;
836 }
837
838 if (pml)
839 {
840 delete omeg_cf;
841 delete omeg2_cf;
842 delete negomeg_cf;
843 delete pml;
844 }
845 delete a;
846 delete q_fec;
847 delete v_fec;
848 delete hatp_fes;
849 delete hatp_fec;
850 delete hatu_fes;
851 delete hatu_fec;
852 delete u_fec;
853 delete p_fec;
854 delete u_fes;
855 delete p_fes;
856
857 return 0;
858}
859
861{
862 return acoustics_solution(x).real();
863}
864
866{
867 return acoustics_solution(x).imag();
868}
869
871{
872 return p_exact_r(X);
873}
874
876{
877 return p_exact_i(X);
878}
879
880void gradp_exact_r(const Vector &x, Vector &grad_r)
881{
882 grad_r.SetSize(x.Size());
883 vector<complex<real_t>> grad;
885 for (unsigned i = 0; i < grad.size(); i++)
886 {
887 grad_r[i] = grad[i].real();
888 }
889}
890
891void gradp_exact_i(const Vector &x, Vector &grad_i)
892{
893 grad_i.SetSize(x.Size());
894 vector<complex<real_t>> grad;
896 for (unsigned i = 0; i < grad.size(); i++)
897 {
898 grad_i[i] = grad[i].imag();
899 }
900}
901
903{
904 return acoustics_solution_laplacian(x).real();
905}
906
908{
909 return acoustics_solution_laplacian(x).imag();
910}
911
912// u = - ∇ p / (i ω )
913// = i (∇ p_r + i * ∇ p_i) / ω
914// = - ∇ p_i / ω + i ∇ p_r / ω
915void u_exact_r(const Vector &x, Vector & u)
916{
917 gradp_exact_i(x,u);
918 u *= -1./omega;
919}
920
921void u_exact_i(const Vector &x, Vector & u)
922{
923 gradp_exact_r(x,u);
924 u *= 1./omega;
925}
926
927void hatu_exact_r(const Vector & X, Vector & hatu)
928{
929 u_exact_r(X,hatu);
930}
931void hatu_exact_i(const Vector & X, Vector & hatu)
932{
933 u_exact_i(X,hatu);
934}
935
936// ∇⋅u = i Δ p / ω
937// = i (Δ p_r + i * Δ p_i) / ω
938// = - Δ p_i / ω + i Δ p_r / ω
939
941{
942 return -d2_exact_i(x)/omega;
943}
944
946{
947 return d2_exact_r(x)/omega;
948}
949
950// f = ∇⋅u + i ω p
951// f_r = ∇⋅u_r - ω p_i
953{
954 real_t p = p_exact_i(x);
955 real_t divu = divu_exact_r(x);
956 return divu - omega * p;
957}
958
959// f_i = ∇⋅u_i + ω p_r
961{
962 real_t p = p_exact_r(x);
963 real_t divu = divu_exact_i(x);
964 return divu + omega * p;
965}
966
967complex<real_t> acoustics_solution(const Vector & X)
968{
969 complex<real_t> zi = complex<real_t>(0., 1.);
970 switch (prob)
971 {
973 case plane_wave:
974 {
975 real_t beta = omega/std::sqrt((real_t)X.Size());
976 complex<real_t> alpha = beta * zi * X.Sum();
977 return exp(alpha);
978 }
979 break;
980 case gaussian_beam:
981 case pml_beam_scatter:
982 {
983 real_t rk = omega;
984 real_t degrees = 45;
985 real_t alpha = (180+degrees) * M_PI/180.;
986 real_t sina = sin(alpha);
987 real_t cosa = cos(alpha);
988 // shift the origin
989 real_t shift = 0.1;
990 real_t xprim=X(0) + shift;
991 real_t yprim=X(1) + shift;
992
993 real_t x = xprim*sina - yprim*cosa;
994 real_t y = xprim*cosa + yprim*sina;
995 //wavelength
996 real_t rl = 2.*M_PI/rk;
997
998 // beam waist radius
999 real_t w0 = 0.05;
1000
1001 // function w
1002 real_t fact = rl/M_PI/(w0*w0);
1003 real_t aux = 1. + (fact*y)*(fact*y);
1004
1005 real_t w = w0*sqrt(aux);
1006
1007 real_t phi0 = atan(fact*y);
1008
1009 real_t r = y + 1./y/(fact*fact);
1010
1011 // pressure
1012 complex<real_t> ze = - x*x/(w*w) - zi*rk*y - zi * real_t(M_PI) * x * x/rl/r +
1013 zi*phi0/2_r;
1014 real_t pf = pow(2.0/M_PI/(w*w),0.25);
1015
1016 return pf*exp(ze);
1017 }
1018 break;
1019 case pml_pointsource:
1020 {
1021 real_t x = X(0)-0.5;
1022 real_t y = X(1)-0.5;
1023 real_t r = sqrt(x*x + y*y);
1024 real_t beta = omega * r;
1025 complex<real_t> Ho = real_t(jn(0, beta)) + zi * real_t(yn(0, beta));
1026 return 0.25_r*zi*Ho;
1027 }
1028 break;
1029 default:
1030 MFEM_ABORT("Should be unreachable");
1031 return 1;
1032 break;
1033 }
1034}
1035
1036void acoustics_solution_grad(const Vector & X, vector<complex<real_t>> & dp)
1037{
1038 dp.resize(X.Size());
1039 complex<real_t> zi = complex<real_t>(0., 1.);
1040 // initialize
1041 for (int i = 0; i<X.Size(); i++) { dp[i] = 0.0; }
1042 switch (prob)
1043 {
1045 case plane_wave:
1046 {
1047 real_t beta = omega/std::sqrt((real_t)X.Size());
1048 complex<real_t> alpha = beta * zi * X.Sum();
1049 complex<real_t> p = exp(alpha);
1050 for (int i = 0; i<X.Size(); i++)
1051 {
1052 dp[i] = zi * beta * p;
1053 }
1054 }
1055 break;
1056 case gaussian_beam:
1057 case pml_beam_scatter:
1058 {
1059 real_t rk = omega;
1060 real_t degrees = 45;
1061 real_t alpha = (180+degrees) * M_PI/180.;
1062 real_t sina = sin(alpha);
1063 real_t cosa = cos(alpha);
1064 // shift the origin
1065 real_t shift = 0.1;
1066 real_t xprim=X(0) + shift;
1067 real_t yprim=X(1) + shift;
1068
1069 real_t x = xprim*sina - yprim*cosa;
1070 real_t y = xprim*cosa + yprim*sina;
1071 real_t dxdxprim = sina, dxdyprim = -cosa;
1072 real_t dydxprim = cosa, dydyprim = sina;
1073 //wavelength
1074 real_t rl = 2.*M_PI/rk;
1075
1076 // beam waist radius
1077 real_t w0 = 0.05;
1078
1079 // function w
1080 real_t fact = rl/M_PI/(w0*w0);
1081 real_t aux = 1. + (fact*y)*(fact*y);
1082
1083 real_t w = w0*sqrt(aux);
1084 real_t dwdy = w0*fact*fact*y/sqrt(aux);
1085
1086 real_t phi0 = atan(fact*y);
1087 real_t dphi0dy = cos(phi0)*cos(phi0)*fact;
1088
1089 real_t r = y + 1./y/(fact*fact);
1090 real_t drdy = 1. - 1./(y*y)/(fact*fact);
1091
1092 constexpr real_t r2 = 2.0;
1093 const real_t rPI = M_PI;
1094
1095 // pressure
1096 complex<real_t> ze = - x*x/(w*w) - zi*rk*y - zi * rPI * x * x/rl/r +
1097 zi*phi0/r2;
1098
1099 complex<real_t> zdedx = -r2*x/(w*w) - r2*zi*rPI*x/rl/r;
1100 complex<real_t> zdedy = r2*x*x/(w*w*w)*dwdy - zi*rk + zi*rPI*x*x/rl/
1101 (r*r)*drdy + zi*dphi0dy/r2;
1102
1103 real_t pf = pow(2.0/M_PI/(w*w),0.25);
1104 real_t dpfdy = -pow(2./M_PI/(w*w),-0.75)/M_PI/(w*w*w)*dwdy;
1105
1106 complex<real_t> zp = pf*exp(ze);
1107 complex<real_t> zdpdx = zp*zdedx;
1108 complex<real_t> zdpdy = dpfdy*exp(ze)+zp*zdedy;
1109
1110 dp[0] = (zdpdx*dxdxprim + zdpdy*dydxprim);
1111 dp[1] = (zdpdx*dxdyprim + zdpdy*dydyprim);
1112 }
1113 break;
1114 default:
1115 MFEM_ABORT("Should be unreachable");
1116 break;
1117 }
1118}
1119
1120complex<real_t> acoustics_solution_laplacian(const Vector & X)
1121{
1122 complex<real_t> zi = complex<real_t>(0., 1.);
1123 switch (prob)
1124 {
1126 case plane_wave:
1127 {
1128 real_t beta = omega/std::sqrt((real_t)X.Size());
1129 complex<real_t> alpha = beta * zi * X.Sum();
1130 return dim * beta * beta * exp(alpha);
1131 }
1132 break;
1133 case gaussian_beam:
1134 case pml_beam_scatter:
1135 {
1136 real_t rk = omega;
1137 real_t degrees = 45;
1138 real_t alpha = (180+degrees) * M_PI/180.;
1139 real_t sina = sin(alpha);
1140 real_t cosa = cos(alpha);
1141 // shift the origin
1142 real_t shift = 0.1;
1143 real_t xprim=X(0) + shift;
1144 real_t yprim=X(1) + shift;
1145
1146 real_t x = xprim*sina - yprim*cosa;
1147 real_t y = xprim*cosa + yprim*sina;
1148 real_t dxdxprim = sina, dxdyprim = -cosa;
1149 real_t dydxprim = cosa, dydyprim = sina;
1150 //wavelength
1151 real_t rl = 2.*M_PI/rk;
1152
1153 // beam waist radius
1154 real_t w0 = 0.05;
1155
1156 // function w
1157 real_t fact = rl/M_PI/(w0*w0);
1158 real_t aux = 1. + (fact*y)*(fact*y);
1159
1160 real_t w = w0*sqrt(aux);
1161 real_t dwdy = w0*fact*fact*y/sqrt(aux);
1162 real_t d2wdydy = w0*fact*fact*(1. - (fact*y)*(fact*y)/aux)/sqrt(aux);
1163
1164 real_t phi0 = atan(fact*y);
1165 real_t dphi0dy = cos(phi0)*cos(phi0)*fact;
1166 real_t d2phi0dydy = -2.*cos(phi0)*sin(phi0)*fact*dphi0dy;
1167
1168 real_t r = y + 1./y/(fact*fact);
1169 real_t drdy = 1. - 1./(y*y)/(fact*fact);
1170 real_t d2rdydy = 2./(y*y*y)/(fact*fact);
1171
1172 constexpr real_t r2 = 2.0;
1173 const real_t rPI = M_PI;
1174
1175 // pressure
1176 complex<real_t> ze = - x*x/(w*w) - zi*rk*y - zi * rPI * x * x/rl/r +
1177 zi*phi0/r2;
1178
1179 complex<real_t> zdedx = -r2*x/(w*w) - r2*zi*rPI*x/rl/r;
1180 complex<real_t> zdedy = r2*x*x/(w*w*w)*dwdy - zi*rk + zi*rPI*x*x/rl/
1181 (r*r)*drdy + zi*dphi0dy/r2;
1182 complex<real_t> zd2edxdx = -r2/(w*w) - r2*zi*rPI/rl/r;
1183 complex<real_t> zd2edxdy = 4_r*x/(w*w*w)*dwdy + zi*r2*rPI*x/rl/(r*r)*drdy;
1184 complex<real_t> zd2edydx = zd2edxdy;
1185 complex<real_t> zd2edydy = -6_r*x*x/(w*w*w*w)*dwdy*dwdy
1186 //+ complex<real_t>(2.*x*x/(w*w*w)*d2wdydy)
1187 + real_t(2.*x*x/(w*w*w)*d2wdydy)
1188 - zi * real_t(2.*M_PI*x*x/rl/(r*r*r)*drdy*drdy)
1189 + zi * real_t(M_PI*x*x/rl/(r*r)*d2rdydy) + zi/real_t(2.*d2phi0dydy);
1190
1191 real_t pf = pow(2.0/M_PI/(w*w),0.25);
1192 real_t dpfdy = -pow(2./M_PI/(w*w),-0.75)/M_PI/(w*w*w)*dwdy;
1193 real_t d2pfdydy = -1./M_PI*pow(2./M_PI,-0.75)*(-1.5*pow(w,-2.5)
1194 *dwdy*dwdy + pow(w,-1.5)*d2wdydy);
1195
1196
1197 complex<real_t> zp = pf*exp(ze);
1198 complex<real_t> zdpdx = zp*zdedx;
1199 complex<real_t> zdpdy = dpfdy*exp(ze)+zp*zdedy;
1200 complex<real_t> zd2pdxdx = zdpdx*zdedx + zp*zd2edxdx;
1201 complex<real_t> zd2pdxdy = zdpdy*zdedx + zp*zd2edxdy;
1202 complex<real_t> zd2pdydx = dpfdy*exp(ze)*zdedx + zdpdx*zdedy + zp*zd2edydx;
1203 complex<real_t> zd2pdydy = d2pfdydy*exp(ze) + dpfdy*exp(
1204 ze)*zdedy + zdpdy*zdedy + zp*zd2edydy;
1205
1206
1207 return (zd2pdxdx*dxdxprim + zd2pdydx*dydxprim)*dxdxprim
1208 + (zd2pdxdy*dxdxprim + zd2pdydy*dydxprim)*dydxprim
1209 + (zd2pdxdx*dxdyprim + zd2pdydx*dydyprim)*dxdyprim
1210 + (zd2pdxdy*dxdyprim + zd2pdydy*dydyprim)*dydyprim;
1211 }
1212 break;
1213 default:
1214 MFEM_ABORT("Should be unreachable");
1215 return 1;
1216 break;
1217 }
1218}
1219
1221{
1222 Vector center(dim);
1223 center = 0.5;
1224 real_t r = 0.0;
1225 for (int i = 0; i < dim; ++i)
1226 {
1227 r += pow(x[i] - center[i], 2.);
1228 }
1229 real_t n = 5.0 * omega / M_PI;
1230 real_t coeff = pow(n, 2) / M_PI;
1231 real_t alpha = -pow(n, 2) * r;
1232 return -omega * coeff * exp(alpha)/omega;
1233}
Dynamic 2D array using row-major layout.
Definition array.hpp:372
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
T & Last()
Return the last element in the array.
Definition array.hpp:802
A class to handle Block diagonal preconditioners in a matrix-free implementation.
Operator & GetDiagonalBlock(int iblock)
Return a reference to block i,i.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
A class to handle Block systems in a matrix-free implementation.
void SetBlock(int iRow, int iCol, Operator *op, real_t c=1.0)
Add a block op in the block-entry (iblock, jblock).
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
int NumRowBlocks() const
Return the number of row blocks.
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
Class for setting up a simple Cartesian PML region.
Definition pml.hpp:19
void SetAttributes(Mesh *mesh_, Array< int > *attrNonPML=nullptr, Array< int > *attrPML=nullptr)
Mark element in the PML region.
Definition pml.cpp:70
void SetOmega(real_t omega_)
Definition pml.hpp:64
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Mimic the action of a complex operator using two real operators.
virtual Operator & imag()
virtual Operator & real()
Real or imaginary part accessor methods.
A coefficient that is constant across space and time.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
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
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
A general function coefficient.
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
The Auxiliary-space Divergence Solver in hypre.
Definition hypre.hpp:1922
The Auxiliary-space Maxwell Solver in hypre.
Definition hypre.hpp:1845
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition hypre.cpp:5111
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Abstract class for hypre's solvers and preconditioners.
Definition hypre.hpp:1162
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
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
Derived matrix coefficient that has the value of the parent matrix coefficient where it is active and...
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
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10558
Geometry::Type GetElementGeometry(int i) const
Definition mesh.hpp:1371
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10626
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void SetIntegrationRule(const IntegrationRule &ir)
Prescribe a fixed IntegrationRule to use.
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.
Class representing the parallel weak formulation. (Convenient for DPG Equations)
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
int GetTrueVSize() const override
Return the number of local vector true dofs.
Definition pfespace.hpp:289
Class for parallel grid function.
Definition pgridfunc.hpp:33
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, const Array< int > &attr)
void MakeRef(FiniteElementSpace *f, real_t *v) override
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Class for parallel meshes.
Definition pmesh.hpp:34
Helper class for ParaView visualization data.
void SetHighOrderOutput(bool high_order_output_)
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
void SetDataFormat(VTKFormat fmt)
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
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
Derived coefficient that takes the value of the parent coefficient for the active attributes and is z...
Matrix coefficient defined as a product of a scalar coefficient and a matrix coefficient.
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:922
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
Vector beta
const real_t alpha
Definition ex15.cpp:369
bool exact_known
Definition ex25.cpp:144
prob_type
Definition ex25.cpp:149
int main()
real_t a
Definition lissajous.cpp:41
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)
real_t detJ_r_function(const Vector &x, CartesianPML *pml)
PML stretching functions: See https://doi.org/10.1006/jcph.1994.1159.
Definition pml.cpp:160
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void Jt_J_detJinv_r_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:191
void abs_Jt_J_detJinv_2_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:223
void Jt_J_detJinv_i_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:207
float real_t
Definition config.hpp:43
real_t abs_detJ_2_function(const Vector &x, CartesianPML *pml)
Definition pml.cpp:180
real_t detJ_i_function(const Vector &x, CartesianPML *pml)
Definition pml.cpp:170
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:486
const char vishost[]
real_t p(const Vector &x, real_t t)
real_t hatp_exact_r(const Vector &X)
void hatu_exact_i(const Vector &X, Vector &hatu)
real_t omega
real_t p_exact_r(const Vector &x)
complex< real_t > acoustics_solution(const Vector &X)
void gradp_exact_i(const Vector &x, Vector &gradu)
void gradp_exact_r(const Vector &x, Vector &gradu)
real_t rhs_func_r(const Vector &x)
complex< real_t > acoustics_solution_laplacian(const Vector &X)
real_t divu_exact_i(const Vector &x)
int dim
void u_exact_r(const Vector &x, Vector &u)
real_t divu_exact_r(const Vector &x)
real_t d2_exact_r(const Vector &x)
real_t source_function(const Vector &x)
real_t hatp_exact_i(const Vector &X)
real_t d2_exact_i(const Vector &x)
real_t rhs_func_i(const Vector &x)
void acoustics_solution_grad(const Vector &X, vector< complex< real_t > > &dp)
prob_type prob
prob_type
@ pml_beam_scatter
@ gaussian_beam
@ plane_wave
@ pml_plane_wave_scatter
@ pml_pointsource
@ pml_general
real_t p_exact_i(const Vector &x)
void hatu_exact_r(const Vector &X, Vector &hatu)
void u_exact_i(const Vector &x, Vector &u)
Helper struct to convert a C++ type to an MPI type.