Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Resonance with square scatterer

More
3 years 4 weeks ago - 3 years 3 weeks ago #3616 by ajf367
Edited: I am very new to NGSolve. I'm trying to verify some resonance values that I've solved for through collocation methods on the integral equation form for the 2D Helmholtz equation with constant value scatterer, the scatterer has value 3, so the equation looks like
[tex]\Delta u + 3 \lambda u =0[/tex]
with Sommerfeld radiation conditions.

I'm using the tutorial as well as things I've seen on other posts and have the following code:
Code:
geo = SplineGeometry() Points = [(0,0), (1,0), (1,1), (0,1)] for pnt in Points: geo.AddPoint(*pnt) for pids in [[0,1],[1,2],[2,3], [3,0]]: geo.Append(["line"]+pids,leftdomain=1,rightdomain=2,bc="scat") geo.AddCircle( (0.5,0.5), 1.2, leftdomain=3, bc="outerbdn") geo.AddCircle( (0.5,0.5), 1.0, leftdomain=2, rightdomain=3, bc="innerbdn") geo.SetMaterial(1, "inner") geo.SetMaterial(2, "air") geo.SetMaterial(3, "pmlregion") mesh = Mesh(geo.GenerateMesh(maxh=0.1)) mesh.Curve(3) mesh.SetPML(pml.Radial(rad=1,alpha=1j,origin=(0.5,0.5)), "pmlregion") omega_0 = 1 omega_tilde = 3 domain_values = {'inner': omega_tilde,'air': omega_0, 'pmlregion': omega_0} values_list = [domain_values[mat] for mat in mesh.GetMaterials()] omega = CoefficientFunction(values_list) Draw(omega, mesh, 'piecewise') fes = H1(mesh, order=4, complex=True, dirichlet="dir") u = fes.TrialFunction() v = fes.TestFunction() a = BilinearForm (fes, symmetric=True) a += grad(u)*grad(v)*dx m = BilinearForm (fes, symmetric=True) m += u*v*dx a.Assemble() m.Assemble() u = GridFunction(fes, multidim=50, name='resonances') with TaskManager(): lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(), u.vecs, shift=0) Draw(u) print ("lam: ", lam)


If someone can let me know if this is actually solving for the resonance values with a [0,1]x[0,1] square scatterer with value 3 I would be greatly appreciative.
Last edit: 3 years 3 weeks ago by ajf367. Reason: Worked some code out and I think if might be better now, still unsure.
Time to create page: 0.134 seconds