search

jackhamel16 / Acoustics

Acoustic integral equation solver
1 stars 0 forks source link

investigate why CFIE doesnt agree well with IE #103

Closed jackhamel16 closed 3 years ago

jackhamel16 commented 3 years ago

The only rime there is strong agreement between the CFIE and non-CFIE is when running on a sphere. For open objects, there is zero agreement. For a cylinder (closed) I get an error around 35% from initial testing. Follow up looking into this.

jackhamel16 commented 3 years ago

cylinder_r1_h2_3201.msh:

l2-norm error between soft IE and soft CFIE without ACA: 36.7% l2-norm error between soft IE and soft CFIE with ACA: 36.7% l2-norm error between soft CFIE w/o ACA and soft CFIE w/ ACA: 2.62e-3% l2-norm error between soft IE w/o ACA and soft IE w/ ACA: 0.0228%

jackhamel16 commented 3 years ago

cube_1m.msh:

l2-norm error between soft IE and soft CFIE without ACA: 32.57% l2-norm error between soft IE and soft CFIE with ACA: 32.57% l2-norm error between soft CFIE w/o ACA and soft CFIE w/ ACA: 4.29e-4% l2-norm error between soft IE w/o ACA and soft IE w/ ACA: 6.44e-3%

jackhamel16 commented 3 years ago

sphere_1m_1266.msh:

l2-norm error between soft IE and soft CFIE without ACA: 1.06% l2-norm error between soft IE and soft CFIE with ACA: 1.06% l2-norm error between soft CFIE w/o ACA and soft CFIE w/ ACA: 1.37e-3% l2-norm error between soft IE w/o ACA and soft IE w/ ACA: 0.02%

jackhamel16 commented 3 years ago

for the tests above, lambda=2pi, spherical wave incident with l=m=0

jackhamel16 commented 3 years ago

Report on this issue is here: Week of 09-13-2020 update.docx Week of 09-13-2020 update.pdf

jackhamel16 commented 3 years ago

It doesn't appear that there is a problem with the elements of a mesh. Using mesh_explorer.py, I can see that the surface elements are all triangles and the number of them matches what my code is seeing as well. The other elements are point elements and line elements, but those are not used in my code. GMSH meshes any mesh by first discretizing nodes with points, then the curves connecting those points with line elements, then those line elements are extended across the surfaces connecting the curves using triangle elements. this can go one step further to mesh with tetrahedral elements, but I stop at 2D since solving integral equations only requires discretizing boundaries

jackhamel16 commented 3 years ago

Investigating direction of normal vectors:

jackhamel16 commented 3 years ago

normal_testing.jl.txt this is the script used to investigate the normal directions. remove .txt extension to use

jackhamel16 commented 3 years ago

Now that I have ruled out the direction of the normals. I will try running problems solving the sound soft normal derivative IE and see if that will match the sound soft IE

jackhamel16 commented 3 years ago

Comparing sound soft IE to sound soft normal derivative IE:

Insights: 1) The only situation out of what has been tested that produces a "reasonable" (really I think it should be much lower, but it is workable with) error is a spherical geometry with a spherical incident wave. I test various values of l and m in the spherical wave and the more spatially varying (i.e. increasing l and m) the wave is, the worse the error becomes. l=m=0: error = 0.0102; l,m=5,-3: error = 0.0244; l,m=8,5: error = 0.0405

jackhamel16 commented 3 years ago

Comparing sound soft IE to sound soft normal derivative IE:

Follow-ups: 1) Comparing sound soft normal deriv IE derivation to coded version and same with the spherical derivative normal derivative.

jackhamel16 commented 3 years ago

Comparing both IEs to analytical solutions for l and m not equal to zero using analytical test framework. Holding Lambda constant for each refinement at the moment since that is what is coded up.

1) Run 3 refinement iterations of sound soft IE and sound soft IE normal derivative with l=m=0 sphere_convergence_results_softIE_l0_m0 sound soft normal derivative below sphere_convergence_results_softIE_l5_m-3

1) Run 3 refinement iterations of sound soft IE and sound soft IE normal derivative with l=5 and m=-3 sphere_convergence_results_softIE_l5_m-3 sphere_convergence_results_softIENormalDeriv_l5_m-3

Conlusions:

The error seems to increase for both the soft IE and soft IE normal derivative in similar amounts when making l and m non zero therefore, I think this further proves that the sphere geometry works fine and other geometries are the problem point.

jackhamel16 commented 3 years ago

How does the soft IE compare to the soft IE normal derivative when the incident wave is a planewave?

The error is not substantially better or worse than when a spherical wave is considered. The error between the two is 0.0164. This could mean that the problem does not lie in the definition of the one of the excitations.

jackhamel16 commented 3 years ago

I went back to reconfirm that the normal directions on the side of the cylinder are correct after realizing maybe my last visualization attempt was incorrect. This was done by using the .pos generating script normal_testing.jl and then checking manually the normal directions of elements assigned 0 (their correct normal directions were not determined by the algorithm in the for loop). Note that GMSH during mesh inspection will give you element tags, which do not directly correspond to ele_idxs in the script. normal_testing.zip

jackhamel16 commented 3 years ago

To check the mesh for any holes, I am going to assign all element sources a value of 1 and visually inspect the results in GMSH.

Doing this, I was unable to notice any holes in the mesh.

jackhamel16 commented 3 years ago

My next test was to compute the volume of the cylinder using the mesh elements as outlined in this document eric meeting sep 21 2021.pdf

The volume of the cylinder being modeled is actual_volume=2pi=6.283 and the volume of the mesh volume=6.272. The relative error between these two: (volume - actual_volume) / volume = -0.00179. I believe this error is negligible because the mesh is created such that the nodes of the mesh exist on the actual geometry being modeled, but the edges connecting the nodes will cut off a small portion of the cylinder (what was an arc becomes a straight line between the nodes) thus minorly reducing the volume of the cylinder. I verified that the nodes on the sides of the cylinder do in fact sit on the actual cylinder by computing their distances from the z axis

As a sanity check, I did this same test with a mesh of 10645 elements and the relative error dropped to -0.00053 which is what I would expect for a finer mesh if what I said above about the meshing process is true

jackhamel16 commented 3 years ago

As a follow on to the volume check above, I will check what the same calculation produces for the sphere with 1266 elements that I know works well with the MFIE and CFIE. The computed volume of the sphere is volume = 4.1516 and the actual volume = 4.1888 creating a relative error = -0.008957.

Now using a sphere with 8010 elements. volume = 4.183 and the relative error = -0.001387

jackhamel16 commented 3 years ago

One thing I have noticed is that the results from the soft IE and soft CFIE do share key features such as singular behavior on edges of the geometry. I am going to now test a geometry that does not have any sharp edges. This will be a cylinder with spherical end caps as pictured below. image Running the volume test on this I get a relative error = -0.007516 which means the normals are correctly pointing outward for all three surfaces.

Spherical Incident Wave

Now running the mesh with a spherical incident wave and the results are as follows with the left image being the magnitude of the surface sources for the soft IE and the right being the magnitude of the surface sources for the soft IE normal derivative: image Obviously, these are not equal. The relative l2-norm of the error between these results is 0.34155. cyl_sph_caps_sph_wave.zip

A plot of the soft IE results divided by soft IE normal derivative results is below with the left plot being the real part and the right being the magnitude: image

Lastly, I want to look at a plot of the relative error of each individual source which is shown below: image

Plane Incident Wave

Now running with a plane wave with lambda=2, k_hat=(0,0,1) (total height of object along z is 4). Left is real part of soft IE and right is real part of soft IE normal derivative: image The relative error here is 0.3453. These results make me want to check if a scale factor can be the problem here.

I divided the real parts of both sets of results and plotted it on the cylinder getting this bizarre result. Left is real part of the scale factors and right is the magnitude: image

Lastly, here is a plot of the relative errors of each element plotted on the object.
image Clearly, for some reason, there bottom cap has correct solutions, but not elsewhere on the object.

jackhamel16 commented 3 years ago

What if I make the cylinder portion of the geometry above much smaller so that the entire object much more closely resembles a sphere? Will the error go down?

The geometry is the same as in the previous test above, but the cylinder height is now 0.1 compared to being 2 previously.

Spherical Incident Wave

l=m=0 and lambda=2. The real part of the surface sources from the soft IE (right) and soft IE normal derivative (left) are plotted below: image There is still significant error here. Since l=m=0, the sources should be nearly constant on this near-sphere object which is what can be seen on the soft IE solution. This is similar to what was seen in the last test with the spherical wave: a spike in the sources at the very ends of the caps (in this test, the spike are also present when looking at the magnitudes instead of real parts). The relative error between these is 0.219.

Here is a plot of the relative error of each element.
image Again, showing that the error spike at the top of the spherical caps. I think two reasonable next tests would be:

1) set cylinder height to 0 and see how that mesh performs. Turns out GMSH won't let you create a cylinder with zero height and I don't want to create one with a very small height because the corresponding mesh would have very tiny elements on the cylinder portion. 2) rotate the object so that the cylinder axis is not on the z-axis and see if anything changes.

jackhamel16 commented 3 years ago

Cylinder height 0.05 Tests

As the title suggests, repeating the last test with a geometry with a cylinder half the size since this is probably as small as I can go without adding a lot more elements. Right is soft IE real part and left is soft IE normal derivative real part image The relative error is 0.0921 which is a definite improvement from the previous tests. A plot of the relative errors tells a similar story as in past tests image

Cylinder height 0.05 Rotated Tests

Now I am going to rotate the object so it aligns with the y axis instead of the z axis. Right is soft IE real and right is soft IE normal derivative real image The relative error is almost identical to when aligned with the z axis at 0.0914. This tells me orientation likely has no effect

jackhamel16 commented 3 years ago

Summary So Far

jackhamel16 commented 3 years ago

Sphere centered off origin

In going through the derivations again for the soft IE normal derivative, I noticed how for the sphere, since it is centered at the origin, the normal vector and r vector will be approximately parallel and for all the other tested geometries, this is not true. If a problem is buried in here, it should show up in this test.
Well, this didn't seem to make a difference. Below are the results of this test for a spherical incident wave with l=m=0 and lambda=2. On the left is real part of soft IE and the right is the real part of the soft IE normal derivative. image The relative error of these two results is 0.01328 -- negligibly different than when centered at the origin.

Sphere with r = 2

not going to bother sharing the results as the rel error is 0.01094. lambda = 3.7 to avoid resonances

Sphere with r = 2 off origin

Combing the last two tests. rel error = 0.01739

Half sphere

Try removing half of a sphere and see if there are problems with that. Ran with a spherical incident wave with l=m=0 and lambda=2.2. Left is real part of soft IE and right is real part of soft IE normal derivative with the top and bottom just showing different angles: image image

A plot of the relative errors of each element from two different angles: image The l2-norm of the relative errors is 0.6583 so things are not agreeing.

What is interesting about this test is that it brings into question the mesh being the source of the problem. The geometry and mesh are created using exactly the same commands with the only difference being that the span of the sphere is set to be Pi instead of 2*Pi so that we only have half of the sphere. Nothing else about the creation of the geometry and mesh is different so GMSH would have to be implicitly changing something based on this one argument to the command that creates the sphere.

Sphere created from scratch

I was originally going to create a union of two half-spheres, but since they technically don't overlap--only share a surface--this didn't work. Instead I created points, a curve connecting the points and extruded the curve to create a sphere. I will test a mesh of this sphere created from scratch and see if the results are any different.

Tested with a spherical incident wave with l=m=0 and lambda=2.2

Finally some interesting results. Left is soft IE real part and right is soft IE normal derivative real part: image The relative error is 3.245. Upon further investigation, it appears I have flipped the normals of the sphere so that they are inward with this construction method. This begs the question of: Why does the normal direction matter here, shouldn't it divide out? Maybe it is because of the principal value of the integral. I don't think it's sign changes with the direction of the normal. This is correct.

I compared the soft IE results with the analytical solution and the rel error is around 0.011 so they are correct.

I inspected the elements in GMSH and can confirm that the nodes are labeled CW when viewed from outside the sphere which will result in inward pointing normals.

Manually flip normals in code

I will rerun the soft IE normal derivative after manually flipping the normal directions in the code. I did this in mesh.jl and adding a -1 times the normals while each one is normalized. Below is the real part of the soft IE normal derivative sources: image While the results are constant, they definitely don't match the soft IE's so there is certainly problem here that only plagues the soft IE normal derivative. Besides the normal vectors, what else about the mesh could only be effecting the soft IE normal derivative? EDIT: Turns out I made a mistake in running the last sphere and with the normals flipped it actually agrees well with the soft IE image the relative error is 0.00772

jackhamel16 commented 3 years ago

Test Results

testing_results.zip

test file for the previous runs: test.jl.txt

jackhamel16 commented 3 years ago

What Now?

Well, it's starting to look like the normals have nothing to do with this bug. I think I need to look closely at the theory and make sure everything lines up with it in code properly with emphasis on parts that are different than what is used in the soft IE.

It appears I may have made an error here in greens_functions.jl by passing [x,y,z]-r_test instead of the negative of that. The math shows that I should pass the negative, but that actually breaks all results, even the sphere.

I changed the code so that I individually pass r_test and r_src to scalarGreensFunctionNormalDerivative and within the function I compute R as R=r_test-r_src (and on the side simplified the calculation of this function).

Now even the sphere test is breaking so I need to retest geometries to try and figure out what is going on starting with a sphere of radius=1 centered at the origin.

Sphere r=1 compared to soft IE

l=m=0 and lambda=3. Below are the surface sources for the soft IE (left) real part and soft IE normal derivative (right) real part. Below them is the same except it shows the imaginary parts and below that is magnitude.
image image image Interestingly, the spread of the source values is way tighter for the soft IE normal derivative now. The rel error here is 1.6438 and that is likely almost entirely due to the fact that the imaginary parts are of opposite sign between the two IEs.

Analytical Comparison

I don't expect it to be accurate, but I know the analytical results are 100% correct so this will be a more concrete place to start investigating compared to using in soft IE results. Use sphere with radius = 1 centered at origin.

lambda = 3

l=m=0 and lambda=3 incident spherical wave.

Below on the left are the real, imaginary and magnitude (top to bottom) of analytical results and right is the same for the soft IE normal derivative results image image image The magnitudes of the sources from the soft IE normal derivative are actually relatively close to the analytical magnitudes. The only noticeable issue is that the imaginary components are relatively equal but opposite in sign.

Below is the results of the analytical sources divided by soft IE normal derivative sources on each element. Left is real part and right is imag part image the average scale factor is -0.450 - 0.957im and it is relatively constant across the sphere. lambda3.zip

lambda = 6

l=m=0 and lambda=6.

Below on the left are the real, imaginary and magnitude (top to bottom) of analytical results and right is the same for the soft IE normal derivative results image image image

Below is the results of the analytical sources divided by soft IE normal derivative sources on each element. Left is real part and right is imag part image The scale factor is quite constant around the sphere. Average scale factor is -2.533 - 2.662im. I don't think there is any significance to the fact that the real and imaginary components of the scale factor are close to each other. lambda6.zip

lambda = 9

l=m=0 and lambda=9.

Below on the left are the real, imaginary and magnitude (top to bottom) of analytical results and right is the same for the soft IE normal derivative results image image image

Below is the results of the analytical sources divided by soft IE normal derivative sources on each element. Left is real part and right is imag part image Again, we see a constant scale factor across the sphere, but it changes with each lambda. Avg scale factor is -5.971 - 4.175im lambda9.zip

lambda=12

This last test I will only share the scale factors only. Left is real part and right is imag part of analytical / soft IE normal derivative: image The avg scale factor is -10.782 - 5.650im lambda12.zip

Summary

Sphere results are no longer valid from soft IE normal derivative (MFIE), but they always seem to differ by some complex scale factor from the analytical solution for a spherical incident wave with l=m=0. The scale factor changes if I change lambda. The scale factor might be dependent on other factors such as radius of the sphere, l, m, etc., but nothing else has been looked at yet.

A plot of the scale factors is shown below: scale_factors_plot The imaginary part actually looks linear while the real part does not. The calculated slope between imag part at lambda=12 and lambda=3 is -0.5214 which when considering I am comparing to the analytical solution with a relatively coarse mesh, this could mean the true scale factors is -0.5.
That being said, I am not sure if that is even helpful right now because the real part does not follow such a linear trend.

Next Steps

I need to try changing other test parameters to determine how this scale factor is effected by the other options. Things to change include, sphere radius, where sphere is centered l and m of spherical wave and swapping spherical wave to plane wave.

jackhamel16 commented 3 years ago

Altering spherical wave l,m

These tests will change l and m of the incident spherical wave with constant lambda =6 on a sphere with radius 1 centered on origin.

l,m=1,0

Below the left spheres are the real, imag, and magnitude (top to bottom) of analytical solution on the sphere. The right spheres are the same but for the soft IE normal derivative solutions. image image image

Below on the left is real part of the analytical sources divided by the soft IE normal derivative sources image The average scale factor is 1.877 + 0.729im

l,m=2,1

Below the left spheres are the real, imag, and magnitude (top to bottom) of analytical solution on the sphere. The right spheres are the same but for the soft IE normal derivative solutions. image image image

Below on the left is real part of the analytical sources divided by the soft IE normal derivative sources image The average scale factor is 1.407 + 0.0514im

Summary

There is a relatively constant scale factor difference between the results of the soft IE normal derivative and analytical solution to the sphere problem. There are a couple elements on the sphere that have a significantly inconsistent scale factor compared to the average. I am not sure why this is as it doesn't seem to be the same elements for each test (for each l,m).

The scale factor is changing as I change l and m. I will move on to another test since l and m are unlikely to be causing the issue, but it is good to have this information in the back of my mind.

jackhamel16 commented 3 years ago

I think I fixed the bug by using the test element's normal vector instead of src element's in computing dG/dn. Need to add tests (left off on fill), update documentation on MFIE. clean everything up. Do verification tests of MFIE. document new files added like in tools.

jackhamel16 commented 3 years ago

Currently, the only tests failing are CFIE ones. I want to verify MFIE is working now before fixing those

Sphere radius 1m

centered at origin with spherical incident wave l,m=2,1 with lambda=6. Convergence test sphere_convergence_results_softIENormalDeriv_l2_m1 Results seem reasonable. Errors converge at rate of-0.9691.

Cylinder with Spherical Caps

Incident plane wave with lambda = 5 traveling along +z. Below the left pills are the real, imag, and magnitude (top to bottom) of analytical solution on the pill. The right pills are the same but for the soft IE normal derivative solutions. image image image The relative error between them is 0.08213. This looks great!

Cube with 1m Edges

Spherical wave l=m=0 lambda=5. Below the left cubes are the real, imag, and magnitude (top to bottom) of analytical solution on the cube. The right cubes are the same but for the soft IE normal derivative solutions. image image image The relative error here is 0.1501 and it looks to me mainly due to the presence of edges. This was ran with a 7 pnt and 1 pnt quadrature rule for src and test integrals. I will try upping the src rule to 13 pnts and see if the error goes down.

Running with 13 points produces a relative error of 0.1523. I would hope this would be converging. 79 point src rule gives an error of 0.1517. Using a 13 pnt source rule and 7 point test rule produces an error of 0.1549. So it appears adding more quadrature points does not force the error to converge.

Cylinder r=1 h=5

spherical incident wave l=m=0, lambda=5. Below the left cylinders are the real, imag, and magnitude (top to bottom) of analytical solution on the cylinder. The right cylinders are the same but for the soft IE normal derivative solutions. image image image The relative error here is worse that it looks sitting at 0.35015. Below is a plot of the magnitude of the relative error of each element image

The error is primarily due to the ends of the cylinders. Overall, the soft IE normal derivative seems to be "smoother" than the soft IE.

Running the same problem with double the number of elements yields a relative error of 0.3522 so there isn't any error convergence occurring as the number of elements increases.

Summary

The bug fixes look to have fixed things for geometries with no edges. There is still a noticeable discrepancy for the pill shape, but I am wondering if that might just be due to the fact that ∂G/∂n→0 when n ̂⋅(r-r^' )→0 which will happen as elements become more and more coplanar. Consider a flat plate. The soft IE normal derivative (MFIE) cannot simulate this structure because n ̂⋅(r-r^' )=0 for all interactions. While this obviously is not the case with the pill, it will have more cancellation than the sphere does and therefore I would not expect the answer to be as close to correct as the soft IE (EFIE). This logic might be flawed, but it is my first thought.

jackhamel16 commented 3 years ago

Bugs fixed. I am now confident the soft IE and soft IE normal derivative are working as they should. Now I need to update the CFIe code

jackhamel16 commented 3 years ago

all tests pass, clean up and merge. follow up checking what old errors were for sphere test vs now.

jackhamel16 commented 3 years ago

This file is full of tests I ran throughout this issue test.jl.txt