MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
acoustics.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11//
12// MFEM Ultraweak DPG example for acoustics (Helmholtz)
13//
14// Compile with: make acoustics
15//
16// Sample runs
17//
18// acoustics -ref 4 -o 1 -rnum 1.0
19// acoustics -m ../../data/inline-tri.mesh -ref 4 -o 2 -sc -rnum 3.0
20// acoustics -m ../../data/amr-quad.mesh -ref 3 -o 3 -sc -rnum 4.5 -prob 1
21// acoustics -m ../../data/inline-quad.mesh -ref 2 -o 4 -sc -rnum 11.5 -prob 1
22// acoustics -m ../../data/inline-hex.mesh -ref 1 -o 2 -sc -rnum 1.0
23
24// Description:
25// This example code demonstrates the use of MFEM to define and solve
26// the "ultraweak" (UW) DPG formulation for the Helmholtz problem
27
28// - Δ p - ω² p = f̃ , in Ω
29// p = p₀, on ∂Ω
30
31// It solves two kinds of problems
32// a) f̃ = 0 and p₀ is a plane wave
33// b) A manufactured solution problem where p_exact is a gaussian beam
34// This example computes and prints out convergence rates for the L² error.
35
36// The DPG UW deals with the First Order System
37// ∇ p + i ω u = 0, in Ω
38// ∇⋅u + i ω p = f, in Ω (1)
39// p = p_0, in ∂Ω
40// where f:=f̃/(i ω)
41
42// Ultraweak-DPG is obtained by integration by parts of both equations and the
43// introduction of trace unknowns on the mesh skeleton
44//
45// p ∈ L²(Ω), u ∈ (L²(Ω))ᵈⁱᵐ
46// p̂ ∈ H^1/2(Ω), û ∈ H^-1/2(Ω)
47// -(p, ∇⋅v) + i ω (u , v) + < p̂, v⋅n> = 0, ∀ v ∈ H(div,Ω)
48// -(u , ∇ q) + i ω (p , q) + < û, q > = (f,q) ∀ q ∈ H¹(Ω)
49// p̂ = p₀ on ∂Ω
50
51// Note:
52// p̂ := p, û := u on the mesh skeleton
53
54// For more information see https://doi.org/10.1016/j.camwa.2017.06.044
55
56// -------------------------------------------------------------
57// | | p | u | p̂ | û | RHS |
58// -------------------------------------------------------------
59// | v | -(p, ∇⋅v) | i ω (u,v) | < p̂, v⋅n> | | |
60// | | | | | | |
61// | q | i ω (p,q) |-(u , ∇ q) | | < û,q > | (f,q) |
62
63// where (q,v) ∈ H¹(Ω) × H(div,Ω)
64
65// Here we use the "Adjoint Graph" norm on the test space i.e.,
66// ||(q,v)||²ᵥ = ||A^*(q,v)||² + ||(q,v)||² where A is the
67// acoustics operator defined by (1)
68
69#include "mfem.hpp"
72#include <fstream>
73#include <iostream>
74
75using namespace std;
76using namespace mfem;
77using namespace mfem::common;
78
79complex<real_t> acoustics_solution(const Vector & X);
80void acoustics_solution_grad(const Vector & X,vector<complex<real_t>> &dp);
81complex<real_t> acoustics_solution_laplacian(const Vector & X);
82
83real_t p_exact_r(const Vector &x);
84real_t p_exact_i(const Vector &x);
85void u_exact_r(const Vector &x, Vector & u);
86void u_exact_i(const Vector &x, Vector & u);
87real_t rhs_func_r(const Vector &x);
88real_t rhs_func_i(const Vector &x);
89void gradp_exact_r(const Vector &x, Vector &gradu);
90void gradp_exact_i(const Vector &x, Vector &gradu);
91real_t divu_exact_r(const Vector &x);
92real_t divu_exact_i(const Vector &x);
93real_t d2_exact_r(const Vector &x);
94real_t d2_exact_i(const Vector &x);
95real_t hatp_exact_r(const Vector & X);
96real_t hatp_exact_i(const Vector & X);
97void hatu_exact_r(const Vector & X, Vector & hatu);
98void hatu_exact_i(const Vector & X, Vector & hatu);
99
100int dim;
102
108
110
111int main(int argc, char *argv[])
112{
113 const char *mesh_file = "../../data/inline-quad.mesh";
114 int order = 1;
115 int delta_order = 1;
116 bool visualization = true;
117 real_t rnum=1.0;
118 int ref = 0;
119 bool static_cond = false;
120 int iprob = 0;
121
122 OptionsParser args(argc, argv);
123 args.AddOption(&mesh_file, "-m", "--mesh",
124 "Mesh file to use.");
125 args.AddOption(&order, "-o", "--order",
126 "Finite element order (polynomial degree)");
127 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
128 "--no-visualization",
129 "Enable or disable GLVis visualization.");
130 args.AddOption(&rnum, "-rnum", "--number_of_wavelengths",
131 "Number of wavelengths");
132 args.AddOption(&iprob, "-prob", "--problem", "Problem case"
133 " 0: plane wave, 1: Gaussian beam");
134 args.AddOption(&delta_order, "-do", "--delta_order",
135 "Order enrichment for DPG test space.");
136 args.AddOption(&ref, "-ref", "--refinements",
137 "Number of serial refinements.");
138 args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
139 "--no-static-condensation", "Enable static condensation.");
140 args.Parse();
141 if (!args.Good())
142 {
143 args.PrintUsage(cout);
144 return 1;
145 }
146 args.PrintOptions(cout);
147
148 if (iprob > 1) { iprob = 0; }
149 prob = (prob_type)iprob;
150
151 omega = 2.*M_PI*rnum;
152
153 Mesh mesh(mesh_file, 1, 1);
154 dim = mesh.Dimension();
155 MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example");
156
157 // Define spaces
158 enum TrialSpace
159 {
160 p_space = 0,
161 u_space = 1,
162 hatp_space = 2,
163 hatu_space = 3
164 };
165 enum TestSpace
166 {
167 q_space = 0,
168 v_space = 1
169 };
170
171 // L2 space for p
172 FiniteElementCollection *p_fec = new L2_FECollection(order-1,dim);
173 FiniteElementSpace *p_fes = new FiniteElementSpace(&mesh,p_fec);
174
175 // Vector L2 space for u
176 FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
177 FiniteElementSpace *u_fes = new FiniteElementSpace(&mesh,u_fec, dim);
178
179 // H^1/2 space for p̂
180 FiniteElementCollection * hatp_fec = new H1_Trace_FECollection(order,dim);
181 FiniteElementSpace *hatp_fes = new FiniteElementSpace(&mesh,hatp_fec);
182
183 // H^-1/2 space for û
184 FiniteElementCollection * hatu_fec = new RT_Trace_FECollection(order-1,dim);
185 FiniteElementSpace *hatu_fes = new FiniteElementSpace(&mesh,hatu_fec);
186
187 // testspace fe collections
188 int test_order = order+delta_order;
189 FiniteElementCollection * q_fec = new H1_FECollection(test_order, dim);
190 FiniteElementCollection * v_fec = new RT_FECollection(test_order-1, dim);
191
192 // Coefficients
193 ConstantCoefficient one(1.0);
194 ConstantCoefficient zero(0.0);
195 Vector vec0(dim); vec0 = 0.;
196 VectorConstantCoefficient vzero(vec0);
197 ConstantCoefficient negone(-1.0);
200 ConstantCoefficient negomeg(-omega);
201
204
205 trial_fes.Append(p_fes);
206 trial_fes.Append(u_fes);
207 trial_fes.Append(hatp_fes);
208 trial_fes.Append(hatu_fes);
209
210 test_fec.Append(q_fec);
211 test_fec.Append(v_fec);
212
213 ComplexDPGWeakForm * a = new ComplexDPGWeakForm(trial_fes,test_fec);
214
215 // i ω (p,q)
216 a->AddTrialIntegrator(nullptr,new MixedScalarMassIntegrator(omeg),
217 TrialSpace::p_space,TestSpace::q_space);
218 // -(u , ∇ q)
219 a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(negone)),
220 nullptr,TrialSpace::u_space,TestSpace::q_space);
221 // -(p, ∇⋅v)
222 a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(one),nullptr,
223 TrialSpace::p_space,TestSpace::v_space);
224 // i ω (u,v)
225 a->AddTrialIntegrator(nullptr,
227 TrialSpace::u_space,TestSpace::v_space);
228 // < p̂, v⋅n>
229 a->AddTrialIntegrator(new NormalTraceIntegrator,nullptr,
230 TrialSpace::hatp_space,TestSpace::v_space);
231 // < û,q >
232 a->AddTrialIntegrator(new TraceIntegrator,nullptr,
233 TrialSpace::hatu_space,TestSpace::q_space);
234
235 // test space integrators (Adjoint graph norm)
236 // (∇q,∇δq)
237 a->AddTestIntegrator(new DiffusionIntegrator(one),nullptr,
238 TestSpace::q_space, TestSpace::q_space);
239 // (q,δq)
240 a->AddTestIntegrator(new MassIntegrator(one),nullptr,
241 TestSpace::q_space, TestSpace::q_space);
242 // (∇⋅v,∇⋅δv)
243 a->AddTestIntegrator(new DivDivIntegrator(one),nullptr,
244 TestSpace::v_space, TestSpace::v_space);
245 // (v,δv)
246 a->AddTestIntegrator(new VectorFEMassIntegrator(one),nullptr,
247 TestSpace::v_space, TestSpace::v_space);
248 // -i ω (∇q,δv)
249 a->AddTestIntegrator(nullptr,new MixedVectorGradientIntegrator(negomeg),
250 TestSpace::q_space, TestSpace::v_space);
251 // i ω (v,∇ δq)
252 a->AddTestIntegrator(nullptr,new MixedVectorWeakDivergenceIntegrator(negomeg),
253 TestSpace::v_space, TestSpace::q_space);
254 // ω^2 (v,δv)
255 a->AddTestIntegrator(new VectorFEMassIntegrator(omeg2),nullptr,
256 TestSpace::v_space, TestSpace::v_space);
257 // - i ω (∇⋅v,δq)
258 a->AddTestIntegrator(nullptr,new VectorFEDivergenceIntegrator(negomeg),
259 TestSpace::v_space, TestSpace::q_space);
260 // i ω (q,∇⋅v)
261 a->AddTestIntegrator(nullptr,new MixedScalarWeakGradientIntegrator(negomeg),
262 TestSpace::q_space, TestSpace::v_space);
263 // ω^2 (q,δq)
264 a->AddTestIntegrator(new MassIntegrator(omeg2),nullptr,
265 TestSpace::q_space, TestSpace::q_space);
266
267 // RHS
270 if (prob == prob_type::gaussian_beam)
271 {
272 a->AddDomainLFIntegrator(new DomainLFIntegrator(f_rhs_r),
273 new DomainLFIntegrator(f_rhs_i), TestSpace::q_space);
274 }
275
280
281 socketstream p_out_r;
282 socketstream p_out_i;
283
284 real_t err0 = 0.;
285 int dof0 = 0; // init to suppress gcc warning
286
287 std::cout << "\n Ref |"
288 << " Dofs |"
289 << " ω |"
290 << " L2 Error |"
291 << " Rate |"
292 << " PCG it |" << endl;
293 std::cout << std::string(60,'-')
294 << endl;
295
296 for (int it = 0; it<=ref; it++)
297 {
298 if (static_cond) { a->EnableStaticCondensation(); }
299 a->Assemble();
300
301 Array<int> ess_tdof_list;
302 Array<int> ess_bdr;
303 if (mesh.bdr_attributes.Size())
304 {
305 ess_bdr.SetSize(mesh.bdr_attributes.Max());
306 ess_bdr = 1;
307 hatp_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
308 }
309
310 // shift the ess_tdofs
311 for (int j = 0; j < ess_tdof_list.Size(); j++)
312 {
313 ess_tdof_list[j] += p_fes->GetTrueVSize() + u_fes->GetTrueVSize();
314 }
315
316 Array<int> offsets(5);
317 offsets[0] = 0;
318 offsets[1] = p_fes->GetVSize();
319 offsets[2] = u_fes->GetVSize();
320 offsets[3] = hatp_fes->GetVSize();
321 offsets[4] = hatu_fes->GetVSize();
322 offsets.PartialSum();
323
324 Vector x(2*offsets.Last());
325 x = 0.;
326
327 GridFunction hatp_gf_r(hatp_fes, x, offsets[2]);
328 GridFunction hatp_gf_i(hatp_fes, x, offsets.Last()+ offsets[2]);
329 hatp_gf_r.ProjectBdrCoefficient(hatpex_r, ess_bdr);
330 hatp_gf_i.ProjectBdrCoefficient(hatpex_i, ess_bdr);
331
332 OperatorPtr Ah;
333 Vector X,B;
334 a->FormLinearSystem(ess_tdof_list,x,Ah, X,B);
335
336 // Setup operator and preconditioner
337 ComplexOperator * Ahc = Ah.As<ComplexOperator>();
338
339 BlockMatrix * A_r = dynamic_cast<BlockMatrix *>(&Ahc->real());
340 BlockMatrix * A_i = dynamic_cast<BlockMatrix *>(&Ahc->imag());
341
342 int num_blocks = A_r->NumRowBlocks();
343 Array<int> tdof_offsets(2*num_blocks+1);
344 tdof_offsets[0] = 0;
345 int k = (static_cond) ? 2 : 0;
346 for (int i=0; i<num_blocks; i++)
347 {
348 tdof_offsets[i+1] = trial_fes[i+k]->GetTrueVSize();
349 tdof_offsets[num_blocks+i+1] = trial_fes[i+k]->GetTrueVSize();
350 }
351 tdof_offsets.PartialSum();
352
353 BlockOperator A(tdof_offsets);
354 for (int i = 0; i<num_blocks; i++)
355 {
356 for (int j = 0; j<num_blocks; j++)
357 {
358 A.SetBlock(i,j,&A_r->GetBlock(i,j));
359 A.SetBlock(i,j+num_blocks,&A_i->GetBlock(i,j), -1.0);
360 A.SetBlock(i+num_blocks,j+num_blocks,&A_r->GetBlock(i,j));
361 A.SetBlock(i+num_blocks,j,&A_i->GetBlock(i,j));
362 }
363 }
364
365 BlockDiagonalPreconditioner M(tdof_offsets);
366 M.owns_blocks = 1;
367 for (int i = 0; i<num_blocks; i++)
368 {
369 M.SetDiagonalBlock(i, new GSSmoother((SparseMatrix&)A_r->GetBlock(i,i)));
370 M.SetDiagonalBlock(num_blocks+i, new GSSmoother((SparseMatrix&)A_r->GetBlock(i,
371 i)));
372 }
373
374 CGSolver cg;
375 cg.SetRelTol(1e-10);
376 cg.SetMaxIter(2000);
377 cg.SetPrintLevel(0);
378 cg.SetPreconditioner(M);
379 cg.SetOperator(A);
380 cg.Mult(B, X);
381
382 a->RecoverFEMSolution(X,x);
383
384 GridFunction p_r(p_fes, x, 0);
385 GridFunction p_i(p_fes, x, offsets.Last());
386
387 GridFunction u_r(u_fes, x, offsets[1]);
388 GridFunction u_i(u_fes, x, offsets.Last() + offsets[1]);
389
392
395
396 int dofs = 0;
397 for (int i = 0; i<trial_fes.Size(); i++)
398 {
399 dofs += trial_fes[i]->GetTrueVSize();
400 }
401
402 real_t p_err_r = p_r.ComputeL2Error(p_ex_r);
403 real_t p_err_i = p_i.ComputeL2Error(p_ex_i);
404 real_t u_err_r = u_r.ComputeL2Error(u_ex_r);
405 real_t u_err_i = u_i.ComputeL2Error(u_ex_i);
406
407 real_t L2Error = sqrt(p_err_r*p_err_r + p_err_i*p_err_i
408 +u_err_r*u_err_r + u_err_i*u_err_i);
409
410 real_t rate_err = (it) ? dim*log(err0/L2Error)/log((real_t)dof0/dofs) : 0.0;
411
412 err0 = L2Error;
413 dof0 = dofs;
414
415 std::ios oldState(nullptr);
416 oldState.copyfmt(std::cout);
417 std::cout << std::right << std::setw(5) << it << " | "
418 << std::setw(10) << dof0 << " | "
419 << std::setprecision(1) << std::fixed
420 << std::setw(4) << 2*rnum << " π | "
421 << std::setprecision(3)
422 << std::setw(10) << std::scientific << err0 << " | "
423 << std::setprecision(2)
424 << std::setw(6) << std::fixed << rate_err << " | "
425 << std::setw(6) << std::fixed << cg.GetNumIterations() << " | "
426 << std::endl;
427 std::cout.copyfmt(oldState);
428
429 if (visualization)
430 {
431 const char * keys = (it == 0 && dim == 2) ? "jRcml\n" : nullptr;
432 char vishost[] = "localhost";
433 int visport = 19916;
434 VisualizeField(p_out_r,vishost, visport, p_r,
435 "Numerical presure (real part)", 0, 0, 500, 500, keys);
436 VisualizeField(p_out_i,vishost, visport, p_i,
437 "Numerical presure (imaginary part)", 501, 0, 500, 500, keys);
438 }
439
440 if (it == ref)
441 {
442 break;
443 }
444
445 mesh.UniformRefinement();
446 for (int i =0; i<trial_fes.Size(); i++)
447 {
448 trial_fes[i]->Update(false);
449 }
450 a->Update();
451 }
452
453 delete a;
454 delete q_fec;
455 delete v_fec;
456 delete hatp_fes;
457 delete hatp_fec;
458 delete hatu_fes;
459 delete hatu_fec;
460 delete u_fec;
461 delete p_fec;
462 delete u_fes;
463 delete p_fes;
464
465 return 0;
466}
467
469{
470 return acoustics_solution(x).real();
471}
472
474{
475 return acoustics_solution(x).imag();
476}
477
479{
480 return p_exact_r(X);
481}
482
484{
485 return p_exact_i(X);
486}
487
488void gradp_exact_r(const Vector &x, Vector &grad_r)
489{
490 grad_r.SetSize(x.Size());
491 vector<complex<real_t>> grad;
493 for (unsigned i = 0; i < grad.size(); i++)
494 {
495 grad_r[i] = grad[i].real();
496 }
497}
498
499void gradp_exact_i(const Vector &x, Vector &grad_i)
500{
501 grad_i.SetSize(x.Size());
502 vector<complex<real_t>> grad;
504 for (unsigned i = 0; i < grad.size(); i++)
505 {
506 grad_i[i] = grad[i].imag();
507 }
508}
509
511{
512
513 return acoustics_solution_laplacian(x).real();
514}
515
517{
518 return acoustics_solution_laplacian(x).imag();
519}
520
521// u = - ∇ p / (i ω )
522// = i (∇ p_r + i * ∇ p_i) / ω
523// = - ∇ p_i / ω + i ∇ p_r / ω
524void u_exact_r(const Vector &x, Vector & u)
525{
526 gradp_exact_i(x,u);
527 u *= -1./omega;
528}
529
530void u_exact_i(const Vector &x, Vector & u)
531{
532 gradp_exact_r(x,u);
533 u *= 1./omega;
534}
535
536void hatu_exact_r(const Vector & X, Vector & hatu)
537{
538 u_exact_r(X,hatu);
539}
540void hatu_exact_i(const Vector & X, Vector & hatu)
541{
542 u_exact_i(X,hatu);
543}
544
545// ∇⋅u = i Δ p / ω
546// = i (Δ p_r + i * Δ p_i) / ω
547// = - Δ p_i / ω + i Δ p_r / ω
549{
550 return -d2_exact_i(x)/omega;
551}
552
554{
555 return d2_exact_r(x)/omega;
556}
557
558// f = ∇⋅u + i ω p
559// f_r = ∇⋅u_r - ω p_i
561{
562 real_t p = p_exact_i(x);
563 real_t divu = divu_exact_r(x);
564 return divu - omega * p;
565}
566
567// f_i = ∇⋅u_i + ω p_r
569{
570 real_t p = p_exact_r(x);
571 real_t divu = divu_exact_i(x);
572 return divu + omega * p;
573}
574
575complex<real_t> acoustics_solution(const Vector & X)
576{
577 complex<real_t> zi = complex<real_t>(0., 1.);
578 switch (prob)
579 {
580 case plane_wave:
581 {
582 real_t beta = omega/std::sqrt((real_t)X.Size());
583 complex<real_t> alpha = beta * zi * X.Sum();
584 return exp(alpha);
585 }
586 break;
587 default:
588 {
589 real_t rk = omega;
590 real_t degrees = 45;
591 real_t alpha = (180+degrees) * M_PI/180.;
592 real_t sina = sin(alpha);
593 real_t cosa = cos(alpha);
594 // shift the origin
595 real_t xprim=X(0) + 0.1;
596 real_t yprim=X(1) + 0.1;
597
598 real_t x = xprim*sina - yprim*cosa;
599 real_t y = xprim*cosa + yprim*sina;
600 //wavelength
601 real_t rl = 2.*M_PI/rk;
602 // beam waist radius
603 real_t w0 = 0.05;
604 // function w
605 real_t fact = rl/M_PI/(w0*w0);
606 real_t aux = 1. + (fact*y)*(fact*y);
607 real_t w = w0*sqrt(aux);
608 real_t phi0 = atan(fact*y);
609 real_t r = y + 1./y/(fact*fact);
610
611 // pressure
612 complex<real_t> ze = - x*x/(w*w) - zi*rk*y - zi * real_t(M_PI) * x * x/rl/r +
613 zi*phi0/2_r;
614
615 real_t pf = pow(2.0/M_PI/(w*w),0.25);
616 return pf*exp(ze);
617 }
618 break;
619 }
620}
621
622void acoustics_solution_grad(const Vector & X, vector<complex<real_t>> & dp)
623{
624 dp.resize(X.Size());
625 complex<real_t> zi = complex<real_t>(0., 1.);
626 switch (prob)
627 {
628 case plane_wave:
629 {
630 real_t beta = omega/std::sqrt((real_t)X.Size());
631 complex<real_t> alpha = beta * zi * X.Sum();
632 complex<real_t> p = exp(alpha);
633 for (int i = 0; i<X.Size(); i++)
634 {
635 dp[i] = zi * beta * p;
636 }
637 }
638 break;
639 default:
640 {
641 real_t rk = omega;
642 real_t degrees = 45;
643 real_t alpha = (180+degrees) * M_PI/180.;
644 real_t sina = sin(alpha);
645 real_t cosa = cos(alpha);
646 // shift the origin
647 real_t xprim=X(0) + 0.1;
648 real_t yprim=X(1) + 0.1;
649
650 real_t x = xprim*sina - yprim*cosa;
651 real_t y = xprim*cosa + yprim*sina;
652 real_t dxdxprim = sina, dxdyprim = -cosa;
653 real_t dydxprim = cosa, dydyprim = sina;
654 //wavelength
655 real_t rl = 2.*M_PI/rk;
656
657 // beam waist radius
658 real_t w0 = 0.05;
659
660 // function w
661 real_t fact = rl/M_PI/(w0*w0);
662 real_t aux = 1. + (fact*y)*(fact*y);
663
664 real_t w = w0*sqrt(aux);
665 real_t dwdy = w0*fact*fact*y/sqrt(aux);
666
667 real_t phi0 = atan(fact*y);
668 real_t dphi0dy = cos(phi0)*cos(phi0)*fact;
669
670 real_t r = y + 1./y/(fact*fact);
671 real_t drdy = 1. - 1./(y*y)/(fact*fact);
672
673 constexpr real_t r2 = 2.0;
674 const real_t rPI = M_PI;
675
676 // pressure
677 complex<real_t> ze = - x*x/(w*w) - zi*rk*y - zi * rPI * x * x/rl/r +
678 zi*phi0/r2;
679
680 complex<real_t> zdedx = -r2*x/(w*w) - r2*zi*rPI*x/rl/r;
681 complex<real_t> zdedy = r2*x*x/(w*w*w)*dwdy - zi*rk + zi*rPI*x*x/rl/
682 (r*r)*drdy + zi*dphi0dy/r2;
683
684 real_t pf = pow(2.0/M_PI/(w*w),0.25);
685 real_t dpfdy = -pow(2./M_PI/(w*w),-0.75)/M_PI/(w*w*w)*dwdy;
686
687 complex<real_t> zp = pf*exp(ze);
688 complex<real_t> zdpdx = zp*zdedx;
689 complex<real_t> zdpdy = dpfdy*exp(ze)+zp*zdedy;
690
691 dp[0] = (zdpdx*dxdxprim + zdpdy*dydxprim);
692 dp[1] = (zdpdx*dxdyprim + zdpdy*dydyprim);
693 if (dim == 3) { dp[2] = 0.0; }
694 }
695 break;
696 }
697}
698
699complex<real_t> acoustics_solution_laplacian(const Vector & X)
700{
701 complex<real_t> zi = complex<real_t>(0., 1.);
702 switch (prob)
703 {
704 case plane_wave:
705 {
706 real_t beta = omega/std::sqrt((real_t)X.Size());
707 complex<real_t> alpha = beta * zi * X.Sum();
708 complex<real_t> p = exp(alpha);
709 return dim * beta * beta * p;
710 }
711 break;
712 default:
713 {
714 real_t rk = omega;
715 real_t degrees = 45;
716 real_t alpha = (180+degrees) * M_PI/180.;
717 real_t sina = sin(alpha);
718 real_t cosa = cos(alpha);
719 // shift the origin
720 real_t xprim=X(0) + 0.1;
721 real_t yprim=X(1) + 0.1;
722
723 real_t x = xprim*sina - yprim*cosa;
724 real_t y = xprim*cosa + yprim*sina;
725 real_t dxdxprim = sina, dxdyprim = -cosa;
726 real_t dydxprim = cosa, dydyprim = sina;
727 //wavelength
728 real_t rl = 2.*M_PI/rk;
729
730 // beam waist radius
731 real_t w0 = 0.05;
732
733 // function w
734 real_t fact = rl/M_PI/(w0*w0);
735 real_t aux = 1. + (fact*y)*(fact*y);
736
737 real_t w = w0*sqrt(aux);
738 real_t dwdy = w0*fact*fact*y/sqrt(aux);
739 real_t d2wdydy = w0*fact*fact*(1. - (fact*y)*(fact*y)/aux)/sqrt(aux);
740
741 real_t phi0 = atan(fact*y);
742 real_t dphi0dy = cos(phi0)*cos(phi0)*fact;
743 real_t d2phi0dydy = -2.*cos(phi0)*sin(phi0)*fact*dphi0dy;
744
745 real_t r = y + 1./y/(fact*fact);
746 real_t drdy = 1. - 1./(y*y)/(fact*fact);
747 real_t d2rdydy = 2./(y*y*y)/(fact*fact);
748
749 constexpr real_t r2 = 2.0;
750 const real_t rPI = M_PI;
751
752 // pressure
753 complex<real_t> ze = - x*x/(w*w) - zi*rk*y - zi * rPI * x * x/rl/r +
754 zi*phi0/r2;
755
756 complex<real_t> zdedx = -r2*x/(w*w) - r2*zi*rPI*x/rl/r;
757 complex<real_t> zdedy = r2*x*x/(w*w*w)*dwdy - zi*rk + zi*rPI*x*x/rl/
758 (r*r)*drdy + zi*dphi0dy/r2;
759 complex<real_t> zd2edxdx = -r2/(w*w) - r2*zi*rPI/rl/r;
760 complex<real_t> zd2edxdy = 4_r*x/(w*w*w)*dwdy + r2*zi*rPI*x/rl/(r*r)*drdy;
761 complex<real_t> zd2edydx = zd2edxdy;
762 complex<real_t> zd2edydy = -6_r*x*x/(w*w*w*w)*dwdy*dwdy + r2*x*x/
763 complex<real_t>(w*w*w)*d2wdydy - r2*zi*rPI*x*x/rl/(r*r*r)*drdy*drdy
764 + zi*rPI*x*x/rl/(r*r)*d2rdydy + zi/r2*d2phi0dydy;
765
766 real_t pf = pow(2.0/M_PI/(w*w),0.25);
767 real_t dpfdy = -pow(2./M_PI/(w*w),-0.75)/M_PI/(w*w*w)*dwdy;
768 real_t d2pfdydy = -1./M_PI*pow(2./M_PI,-0.75)*(-1.5*pow(w,-2.5)
769 *dwdy*dwdy + pow(w,-1.5)*d2wdydy);
770
771 complex<real_t> zp = pf*exp(ze);
772 complex<real_t> zdpdx = zp*zdedx;
773 complex<real_t> zdpdy = dpfdy*exp(ze)+zp*zdedy;
774 complex<real_t> zd2pdxdx = zdpdx*zdedx + zp*zd2edxdx;
775 complex<real_t> zd2pdxdy = zdpdy*zdedx + zp*zd2edxdy;
776 complex<real_t> zd2pdydx = dpfdy*exp(ze)*zdedx + zdpdx*zdedy + zp*zd2edydx;
777 complex<real_t> zd2pdydy = d2pfdydy*exp(ze) + dpfdy*exp(
778 ze)*zdedy + zdpdy*zdedy + zp*zd2edydy;
779
780
781 return (zd2pdxdx*dxdxprim + zd2pdydx*dydxprim)*dxdxprim
782 + (zd2pdxdy*dxdxprim + zd2pdydy*dydxprim)*dydxprim
783 + (zd2pdxdx*dxdyprim + zd2pdydx*dydyprim)*dxdyprim
784 + (zd2pdxdy*dxdyprim + zd2pdydy*dydyprim)*dydyprim;
785 }
786 break;
787 }
788}
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 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
@ gaussian_beam
@ plane_wave
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)
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.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
int NumRowBlocks() const
Return the number of row blocks.
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).
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 representing the DPG weak formulation for complex valued systems (see the class DPGWeakForm).
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.
for Raviart-Thomas elements
Class for domain integration .
Definition lininteg.hpp:109
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:588
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
A general function coefficient.
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
void ProjectBdrCoefficient(Coefficient &coeff, const Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Definition gridfunc.hpp:469
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:322
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition solvers.hpp:260
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:445
Data type sparse matrix.
Definition sparsemat.hpp:51
Vector coefficient that is constant in space and time.
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1283
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
Vector beta
const real_t alpha
Definition ex15.cpp:369
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)
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
const char vishost[]
real_t p(const Vector &x, real_t t)