jackhamel16 / acoustics Goto Github PK

This function doesn't need to be part of the main code. Store it in a new directory called 'tools' and it can be ran when needed and ask for a file name containing WS time delays.

use multiple dispatch to write another version that compute scattering matrix and its derivative in the function rather than take them as arguments.

inside the buildPulseMesh function, there is an indexing into elements during the reshaping done here which is where the error occurs.

ERROR: BoundsError: attempt to access 0-element Array{Array{UInt64,1},1} at index [2]

This is happening because gmsh module is not properly reading the .msh file. If you look at both nodes and elements, they contain no information. This was done on windows. test.jl in packages is set up for testing this problem.

look for any places to speed it up. maybe changing Vs_trans to be Vs thus becoming column major when indexed.

Theory is complete and can now try implementing this

I need the ability to output my sources to a file format friendly to Gmsh for visualization. Use Utkarsh's code to get started. It will probably require some modification.

add a function to compute the interactions between elements in two nodes (provided as inputs to the function) using ACA for the sound-soft IE.

1. Figure out what quadrature rules are of interest right now
2. Build generalized quadrature routine
3. Add rules

This object will have coordinates of vertices and centroid. Use consistent ordering of vertices. Go clockwise.

I need to group elements in sub-blocks (cubes). ACA requires knowledge of well-separated groups of elements to compress the sub Z matrix for those groups to accelerate calculation of the sources. This will be done as follows

Test 1
create a single box encapsulating all elements of a simple geometry (maybe a square plate with 2 elements centered at origin in XY-plane). Check that all elements are within this box.

Test 2
For the same geometry, sort into 8 cubes that are subdivisions of the single cube. Only two should be full and only contain one element each. This will be an edge case since the geometry will lie exactly on the interface of two boxes.

Accomplishing this will first require a data structure to store the octree. This is simplified by the fact that I am first implementing this with the naïve approach meaning I am not going to even use an octree yet. I will partition the space into cubes at my desired granularity, then loop through all elements to determine which cube each lies in and store the index of the element in an array associated with each cube. Once these tests pass, then refactor to use a recursive octree construction algorithm. Details on how to do so will follow later.

Two routines to add

1. Filling the impedance matrix
2. Filling the RHS vector

For the incident wave, use a plane wave
for matrix and RHS do galerkin testing with 1 point rule (point matching) or galerkin (more than one point)

This is for the first iteration of the main routine and will do the following.

• read in mesh
• define incident field as a plane wave
• fill matrix and rhs
• invert to solve for sources

Outputting data will be a follow on issue

The Compression Ratio (CR) is defined as (sum of all length(sub-Z matrices)) / length(direct Z) = CR. If CR > 1, then ACA should not be used to approximate the sub-Z matrix between those nodes. I think this is best added inside the if-else above here where the size of U and V, if ACA is used, or the size of the direct sub-Z is added to a running total outside the for loops tracking the total number of entries representing the compressed Z matrix.

Record rank of each compressed matrix and compute mean, min and max of the ranks

Record number of elements in each node and compute the mean, min and max of the number of elements per node.

On second thought, this could all be done in post-processing in a separate function as follows:

create array for number of elements in nodes size=number of nodes.
create empty array for ranks
initialize variable to count number of matrix elements stored
loop through each node
store number of elements in node
if compressed:
add size(U) + size(V) to counter variable
store rank which is number of columns of U or rows of V
else:
add size(Z) to counter variable
compute the means, mins and maxes as stated above

Usage of static arrays is not clear yet, but was recommended by Eric. Parameters package could potentially clean up the passing of mesh parameters to functions.

run sphere and plate geometries with different setting to determine what settings work best for what types of problems.

Z matrix is used repeatedly in post-processing and time would be saved by storing it

this will include only plane waves for the moment

This will take a sphere centered at the origin with arbitrary radius and compute the pressure sources on the sphere for an arbitrary incident spherical wave.

Just run the test with all four meshes overnight and update last reference convergence rate (slope).

in solving at WS modes and in contructing the scattering matrix

It seems like files are being included multiple times (probably once for each file they're called in). Look into preventing unnecessary includes.

https://discourse.julialang.org/t/does-julia-have-something-like-ifdef-in-c-for-conditional-compile/44721

arrays are row major which impacts performance when indexing them and that is done throughout the code.

The quadrature points on each triangle are computed with every call to integrateTriangle. Many calls to integrateTriangle are integrating over the same triangle at the same points so I am repeating this calculation redundantly. By moving it to buildPulseMesh and storing the quadrature points as members of PulseMesh types I should be able to save computation time.

It also might be smart to add in the parameters package during this refactor.

Another important test for the CFIE is to see if it can resolve resonances over a frequency range.

1. Determine resonances for a sphere analytically
2. Run CFIE, MFIE and EFIE codes sweeping over one or more resonance frequencies
3. Inspect sources from each code.
4. Plot condition number of Z for each IE and each frequency. Condition number should be lowest for CFIE near resonances.
5. option: Plot determinant of Z. If zero, there can be nonzero sources for zero incident fields (modes exist in cavity)

According to my limited knowledge of the profiler and my understanding of my code, the most time is probably spent in the matrix fill routine and sub routines of it such as quadrature. Try out some speed ups on the matrix fill routine and see what happens!

The sound soft IE is only required for construction of the scattering matrix and WS matrix. Once that is done, and the desired excitation for the WS mode is constructed, any of the IE can be solved with that excitation. I could go in and add this feature to my code, if it becomes necessary.

All of these functions need to be tested for type issues, optimized for performance and any general refactoring deemed necessary. Also make sure they all have headings describing the functions' purposes.

create rectangular plate with 8 elements that has four axes of symmetry. Test code to make sure symmetry is in results.

For singular integration routine, do 1/R, over entire triangle, not sub-triangles.

Add a function that will parse a text file containing all of the desired settings for running an instance of the acoustic code.

It seems like there are useless members of the PulseMesh struct. I think only nodes, elements and num_elements are used for anything, but this needs to be verified before refactoring.

Add ability to integrate singularities

This is now done twice so it should be moved to somewhere independent of the functions that need this information and passed as an argument to them.

This file will be called from the command line and take a command line argument specifying the input file name, call the parsing function on that file created in issue #81, then run the proper solver using input settings.

1. Compute Q for geometry of interest making sure to use enough harmonics for scattering matrix
2. Diagonalize Q. Need to test to make sure the diagonalization matches the theory
3. Create new incident field from first column of W. Probably will write a new incident field function that sums up a bunch of incident fields weighted by the eigenvector entries
4. Solve sound soft IE using the new WS mode incident field using the original methods since Z is not stored

When reading in a mesh, if the nodes were not read in in numerical order based on tags, then the elements would not get the correct node locations that comprise themselves. This was discovered when looking over Utkarsh's code for post processing and comparing it to my preprocessing code.

Need to change preprocessing to account for unordered nodes and place them in their correct places in the node locations array that elements uses.

Add post-processing routine to compute scattering matrix as detailed in documentation.

When do you truncate the harmonics used?
Max value of l should be proportional to ka where a is the radius of the scatterer. Some tests to see this are:

1. Plot Bessel function of order l with argument ka. as l increases for a given ka, it should flatten (i.e. the V for that l will become negligible)
2. Plot l2-norm of V vector as l goes up and it should go to zero as l increases. When it is sufficiently small, the sum can be truncated at that l value.

Testing

First test should check if the scattering matrix is symmetric (it should be). Second test is if S is such that S \dot adjoint(S) = 1

Follow derivations in write up.

Test by approximating derivative with finite difference

Currently, a function to evaluate over the integration domain is passed to the gaussQuadrature function as an argument. This is akin to something that would be done in C++, but may not be optimal for Julia. Think of another way to do this and benchmark as it may result in cleaner code and/or a performance boost.

Maybe it could be as simple as evaluating the function at the quadrature points outside of the function then passing the results in as an array?

Put everything together so that ACA is used to compute Z for well separated elements and Z is computed normally for close elements. Use an iterative solver to find the unknowns, J.

1. set up mesh and octree structs
2. for each test node
   for each src node
   if nodes farther than some tolerance, use ACA to compute decomposed sub-Z, U and V
   else nodes are too close, use direct calculations of sub-Z
3. call iterative solver on function that computes Z*J implicitly

Number 2 should be wrapped in a routine that does an effective matrix vector product between Z and and given J.

Question: Should I be storing all these sub-Z matrices (compressed and uncompressed) so that I am not recomputing them with each iteration of the iterative solve? Probably...

Maybe store with each node an array of the sub-Z matrix for that node's interaction with every other node. For uncompressed, it will be a single Z matrix, for compressed it will be an array of two arrays–one U and one V. Now I can write a function to compute ZJ utilizing multiple dispatch. If the provided Z is one matrix, then it will do ZJ, if the provided Z is an array of two matrices, it will compute U*(V*J)

Split this up into multiple tasks:

First look into iterative solvers
1: add ability to store Z and U,V for each node interaction
2: add routine to compute matvec using U and V or Z
3: wrap matvec routine so that given a vector J, it loops through each node pair and computes the corresponding contributions to V and adds them up to produce the entire V vector.

Lastly, follow original steps 1-3 where step three uses the new matvec function

I need an input deck to enable running my code.

Use Utkarsh's routine as a starting point. Need to write tests for this. Try starting with the textbook to generate reference solutions.

should first check if singular and then not compute distance from the source triangle since that is only necessary if it might be a near source integration.

given a vectors J and V and two nodes of an octree, compute the contributions to V in ZJ=V for the elements in the nodes for the sound-soft integral equation.

An input file to specify input parameters and equations to use needs to be created and an output file with source's and their locations. Also need to add options to compute scattering matrix.

computing the triangle area is done for each matrix and vector entry in the fill routines. It should be done once for each triangle and stored in the PulseMesh object.