diff --git a/.travis.yml b/.travis.yml index 12296e7..c9964a4 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,23 +4,26 @@ os: - linux - osx julia: - - 0.6 + - 0.7 branches: - except: - - dev + only: + - master + - /^v\d+\.\d+(\.\d+)?(-\S*)?$/ notifications: email: false env: - PYTHON=Conda-python -script: - - export LD_LIBRARY_PATH=$HOME/.julia/v0.6/Conda/deps/usr/lib; LD_PRELOAD=${HOME}/.julia/v0.6/Conda/deps/usr/lib/libz.so julia -e 'using PyCall;pygui(:tk);using PyPlot;Pkg.clone(pwd());Pkg.test("ParticleScattering"; coverage=true);' -before_install: - - julia -e 'ENV["PYTHON"]=""; Pkg.add("Conda"); using Conda; Conda.update(); Conda.add("matplotlib"); Conda.add("basemap"); Pkg.add("PyCall"); Pkg.build("PyCall"); Pkg.add("PyPlot");' after_success: - #doc - - julia -e 'Pkg.add("Documenter")' - - julia -e 'cd(Pkg.dir("ParticleScattering")); include(joinpath("docs", "make.jl"))' - # push coverage results to Coveralls - # - julia -e 'cd(Pkg.dir("ParticleScattering")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())' # push coverage results to Codecov - - julia -e 'cd(Pkg.dir("ParticleScattering")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())' + - julia -e 'Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())' + +jobs: + include: + - stage: "Documentation" + julia: 0.7 + os: linux + script: + - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); + Pkg.instantiate()' + - julia --project=docs/ docs/make.jl + after_success: skip diff --git a/Manifest.toml b/Manifest.toml new file mode 100644 index 0000000..c98c3c2 --- /dev/null +++ b/Manifest.toml @@ -0,0 +1,471 @@ +[[Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[BinDeps]] +deps = ["Compat", "Libdl", "SHA", "URIParser"] +git-tree-sha1 = "12093ca6cdd0ee547c39b1870e0c9c3f154d9ca9" +uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee" +version = "0.8.10" + +[[BinaryProvider]] +deps = ["Libdl", "SHA"] +git-tree-sha1 = "c7361ce8a2129f20b0e05a89f7070820cfed6648" +uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" +version = "0.5.4" + +[[Blosc]] +deps = ["BinaryProvider", "CMakeWrapper", "Compat", "Libdl"] +git-tree-sha1 = "71fb23581e1f0b0ae7be8ccf0ebfb3600e23ca41" +uuid = "a74b3585-a348-5f62-a45c-50e91977d574" +version = "0.5.1" + +[[BufferedStreams]] +deps = ["Compat", "Test"] +git-tree-sha1 = "5d55b9486590fdda5905c275bb21ce1f0754020f" +uuid = "e1450e63-4bb3-523b-b2a4-4ffa8c0fd77d" +version = "1.0.0" + +[[CMake]] +deps = ["BinDeps", "Libdl", "Test"] +git-tree-sha1 = "6e39bef3cbb8321e8a464b18a5c20d7cef813938" +uuid = "631607c0-34d2-5d66-819e-eb0f9aa2061a" +version = "1.1.1" + +[[CMakeWrapper]] +deps = ["BinDeps", "CMake", "Libdl", "Parameters", "Test"] +git-tree-sha1 = "16d4acb3d37dc05b714977ffefa8890843dc8985" +uuid = "d5fb7624-851a-54ee-a528-d3f3bac0b4a0" +version = "0.2.3" + +[[CSV]] +deps = ["CategoricalArrays", "DataFrames", "DataStreams", "Dates", "Mmap", "Parsers", "Profile", "Random", "Tables", "Test", "Unicode", "WeakRefStrings"] +git-tree-sha1 = "b92c6f626a044cc9619156d54994b94084d40abe" +uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +version = "0.4.3" + +[[Calculus]] +deps = ["Compat"] +git-tree-sha1 = "f60954495a7afcee4136f78d1d60350abd37a409" +uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" +version = "0.4.1" + +[[CategoricalArrays]] +deps = ["Compat", "Future", "Missings", "Printf", "Reexport", "Requires"] +git-tree-sha1 = "94d16e77dfacc59f6d6c1361866906dbb65b6f6b" +uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597" +version = "0.5.2" + +[[CodecZlib]] +deps = ["BinaryProvider", "Libdl", "Test", "TranscodingStreams"] +git-tree-sha1 = "36bbf5374c661054d41410dc53ff752972583b9b" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.5.2" + +[[ColorTypes]] +deps = ["FixedPointNumbers", "Random"] +git-tree-sha1 = "10050a24b09e8e41b951e9976b109871ce98d965" +uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +version = "0.8.0" + +[[Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Printf", "Reexport", "Test"] +git-tree-sha1 = "9f0a0210450acb91c730b730a994f8eef1d3d543" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.9.5" + +[[CommonSubexpressions]] +deps = ["Test"] +git-tree-sha1 = "efdaf19ab11c7889334ca247ff4c9f7c322817b0" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.2.0" + +[[Compat]] +deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] +git-tree-sha1 = "84aa74986c5b9b898b0d1acaf3258741ee64754f" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "2.1.0" + +[[Conda]] +deps = ["Compat", "JSON", "VersionParsing"] +git-tree-sha1 = "b625d802587c2150c279a40a646fba63f9bd8187" +uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d" +version = "1.2.0" + +[[DataFrames]] +deps = ["CategoricalArrays", "CodecZlib", "Compat", "DataStreams", "Dates", "InteractiveUtils", "IteratorInterfaceExtensions", "LinearAlgebra", "Missings", "Printf", "Random", "Reexport", "SortingAlgorithms", "Statistics", "StatsBase", "TableTraits", "Tables", "Test", "TranscodingStreams", "Unicode", "WeakRefStrings"] +git-tree-sha1 = "9cfed75401d25d281076eb5d82de148ac2933f9e" +uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +version = "0.17.1" + +[[DataStreams]] +deps = ["Dates", "Missings", "Test", "WeakRefStrings"] +git-tree-sha1 = "69c72a1beb4fc79490c361635664e13c8e4a9548" +uuid = "9a8bc11e-79be-5b39-94d7-1ccc349a1a85" +version = "0.4.1" + +[[DataStructures]] +deps = ["InteractiveUtils", "OrderedCollections", "Random", "Serialization", "Test"] +git-tree-sha1 = "ca971f03e146cf144a9e2f2ce59674f5bf0e8038" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.15.0" + +[[Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[DelimitedFiles]] +deps = ["Mmap"] +uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" + +[[DiffEqDiffTools]] +deps = ["LinearAlgebra", "Test"] +git-tree-sha1 = "30f82c63cb9d293513b360572e283c19577fcf82" +uuid = "01453d9d-ee7c-5054-8395-0335cb756afa" +version = "0.8.1" + +[[DiffResults]] +deps = ["Compat", "StaticArrays"] +git-tree-sha1 = "34a4a1e8be7bc99bc9c611b895b5baf37a80584c" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "0.0.4" + +[[DiffRules]] +deps = ["Random", "Test"] +git-tree-sha1 = "dc0869fb2f5b23466b32ea799bd82c76480167f7" +uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" +version = "0.0.10" + +[[Distributed]] +deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[FileIO]] +deps = ["Pkg", "Random", "Test"] +git-tree-sha1 = "da32159d4a2e526338506685e280e39ed2f18961" +uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +version = "1.0.6" + +[[FixedPointNumbers]] +deps = ["Test"] +git-tree-sha1 = "b8045033701c3b10bf2324d7203404be7aef88ba" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.5.3" + +[[ForwardDiff]] +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"] +git-tree-sha1 = "4c4d727f1b7e0092134fabfab6396b8945c1ea5b" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "0.10.3" + +[[Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" + +[[HDF5]] +deps = ["BinDeps", "Blosc", "Homebrew", "Libdl", "Mmap", "WinRPM"] +git-tree-sha1 = "e6f0c154d01faef0d0831d075aa8f279f95946da" +uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +version = "0.11.1" + +[[HTTPClient]] +deps = ["Compat", "LibCURL"] +git-tree-sha1 = "161d5776ae8e585ac0b8c20fb81f17ab755b3671" +uuid = "0862f596-cf2d-50af-8ef4-f2be67dfa83f" +version = "0.2.1" + +[[Homebrew]] +deps = ["BinDeps", "InteractiveUtils", "JSON", "Libdl", "Test", "Unicode"] +git-tree-sha1 = "f01fb2f34675f9839d55ba7238bab63ebd2e531e" +uuid = "d9be37ee-ecc9-5288-90f1-b9ca67657a75" +version = "0.7.1" + +[[InteractiveUtils]] +deps = ["LinearAlgebra", "Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[IterativeSolvers]] +deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays", "Test"] +git-tree-sha1 = "a7b8056ecd43aa0d112d4b5d04180194d12c3c62" +uuid = "42fd0dbc-a981-5370-80f2-aaf504508153" +version = "0.7.1" + +[[IteratorInterfaceExtensions]] +git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" +uuid = "82899510-4779-5014-852e-03e436cf321d" +version = "1.0.0" + +[[JLD]] +deps = ["Compat", "FileIO", "HDF5", "LegacyStrings", "Profile", "Random"] +git-tree-sha1 = "95fd5d7f129918a75d0535aaaf5b8e235e6e0b0b" +uuid = "4138dd39-2aa7-5051-a626-17a0bb65d9c8" +version = "0.9.1" + +[[JSON]] +deps = ["Dates", "Distributed", "Mmap", "Sockets", "Test", "Unicode"] +git-tree-sha1 = "1f7a25b53ec67f5e9422f1f551ee216503f4a0fa" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.20.0" + +[[LaTeXStrings]] +deps = ["Compat"] +git-tree-sha1 = "7ab9b8788cfab2bdde22adf9004bda7ad9954b6c" +uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +version = "1.0.3" + +[[LegacyStrings]] +deps = ["Compat"] +git-tree-sha1 = "d4b9bde2694c552fe579cc4462733f1ce08733fe" +uuid = "1b4a561d-cfcb-5daf-8433-43fcf8b4bea3" +version = "0.4.1" + +[[LibCURL]] +deps = ["BinaryProvider", "Libdl"] +git-tree-sha1 = "5ee138c679fa202ebe211b2683d1eee2a87b3dbe" +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.5.1" + +[[LibExpat]] +deps = ["Compat"] +git-tree-sha1 = "fde352ec13479e2f90e57939da2440fb78c5e388" +uuid = "522f3ed2-3f36-55e3-b6df-e94fee9b0c07" +version = "0.5.0" + +[[LibGit2]] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[Libz]] +deps = ["BufferedStreams", "Random", "Test"] +git-tree-sha1 = "d405194ffc0293c3519d4f7251ce51baac9cc871" +uuid = "2ec943e9-cfe8-584d-b93d-64dcb6d567b7" +version = "1.0.0" + +[[LineSearches]] +deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf", "Test"] +git-tree-sha1 = "54eb90e8dbe745d617c78dee1d6ae95c7f6f5779" +uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +version = "7.0.1" + +[[LinearAlgebra]] +deps = ["Libdl"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[LinearMaps]] +deps = ["LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "ea9bb3d793917bced9c98da5ed8b6f449c4a74ba" +uuid = "7a12625a-238d-50fd-b39a-03d52299707e" +version = "2.3.0" + +[[Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[MacroTools]] +deps = ["Compat"] +git-tree-sha1 = "3fd1a3022952128935b449c33552eb65895380c1" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.4.5" + +[[Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[Missings]] +deps = ["Dates", "InteractiveUtils", "SparseArrays", "Test"] +git-tree-sha1 = "d1d2585677f2bd93a97cfeb8faa7a0de0f982042" +uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +version = "0.4.0" + +[[Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[NLSolversBase]] +deps = ["Calculus", "DiffEqDiffTools", "DiffResults", "Distributed", "ForwardDiff", "LinearAlgebra", "Random", "SparseArrays", "Test"] +git-tree-sha1 = "0c6f0e7f2178f78239cfb75310359eed10f2cacb" +uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" +version = "7.3.1" + +[[NaNMath]] +deps = ["Compat"] +git-tree-sha1 = "ce3b85e484a5d4c71dd5316215069311135fa9f2" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "0.3.2" + +[[Optim]] +deps = ["Calculus", "Compat", "DiffEqDiffTools", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "Random", "SparseArrays", "StatsBase", "Test"] +git-tree-sha1 = "0ed8809e72d3b595930c8ade539838a99f94dc8b" +uuid = "429524aa-4258-5aef-a3af-852621145aeb" +version = "0.16.0" + +[[OrderedCollections]] +deps = ["Random", "Serialization", "Test"] +git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.1.0" + +[[Parameters]] +deps = ["Markdown", "OrderedCollections", "REPL", "Test"] +git-tree-sha1 = "70bdbfb2bceabb15345c0b54be4544813b3444e4" +uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" +version = "0.10.3" + +[[Parsers]] +deps = ["Dates", "Mmap", "Test", "WeakRefStrings"] +git-tree-sha1 = "d4e634c9cab597f0ec8f5a1d62c0787e98602c55" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "0.2.22" + +[[Pkg]] +deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" + +[[PositiveFactorizations]] +deps = ["LinearAlgebra", "Test"] +git-tree-sha1 = "957c3dd7c33895469760ce873082fbb6b3620641" +uuid = "85a6dd25-e78a-55b7-8502-1745935b8125" +version = "0.2.2" + +[[Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[Profile]] +deps = ["Printf"] +uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" + +[[PyCall]] +deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Pkg", "Serialization", "Statistics", "Test", "VersionParsing"] +git-tree-sha1 = "6e5bac1b1faf3575731a6a5b76f638f2389561d3" +uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +version = "1.91.2" + +[[PyPlot]] +deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"] +git-tree-sha1 = "3ea391c3211f799df20396d256cee5c8bd755605" +uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee" +version = "2.8.1" + +[[REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[Random]] +deps = ["Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[RecipesBase]] +deps = ["Random", "Test"] +git-tree-sha1 = "b18c3d875ad6805e3db59c4e81f206f04df71d66" +uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +version = "0.5.0" + +[[Reexport]] +deps = ["Pkg"] +git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "0.2.0" + +[[Requires]] +deps = ["Test"] +git-tree-sha1 = "f6fbf4ba64d295e146e49e021207993b6b48c7d1" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "0.5.2" + +[[SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" + +[[Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[SortingAlgorithms]] +deps = ["DataStructures", "Random", "Test"] +git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd" +uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" +version = "0.3.1" + +[[SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[SpecialFunctions]] +deps = ["BinDeps", "BinaryProvider", "Libdl", "Test"] +git-tree-sha1 = "0b45dc2e45ed77f445617b99ff2adf0f5b0f23ea" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "0.7.2" + +[[StaticArrays]] +deps = ["InteractiveUtils", "LinearAlgebra", "Random", "Statistics", "Test"] +git-tree-sha1 = "3841b39ed5f047db1162627bf5f80a9cd3e39ae2" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "0.10.3" + +[[Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[StatsBase]] +deps = ["DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"] +git-tree-sha1 = "8a0f4b09c7426478ab677245ab2b0b68552143c7" +uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +version = "0.30.0" + +[[TableTraits]] +deps = ["IteratorInterfaceExtensions"] +git-tree-sha1 = "b1ad568ba658d8cbb3b892ed5380a6f3e781a81e" +uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" +version = "1.0.0" + +[[Tables]] +deps = ["IteratorInterfaceExtensions", "LinearAlgebra", "Requires", "TableTraits", "Test"] +git-tree-sha1 = "719d5be11e89ae29d79b469c4238b63b53544d38" +uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +version = "0.1.18" + +[[Test]] +deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[TranscodingStreams]] +deps = ["Random", "Test"] +git-tree-sha1 = "a25d8e5a28c3b1b06d3859f30757d43106791919" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.9.4" + +[[URIParser]] +deps = ["Test", "Unicode"] +git-tree-sha1 = "6ddf8244220dfda2f17539fa8c9de20d6c575b69" +uuid = "30578b45-9adc-5946-b283-645ec420af67" +version = "0.4.0" + +[[UUIDs]] +deps = ["Random"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[VersionParsing]] +deps = ["Compat"] +git-tree-sha1 = "c9d5aa108588b978bd859554660c8a5c4f2f7669" +uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" +version = "1.1.3" + +[[WeakRefStrings]] +deps = ["Missings", "Random", "Test"] +git-tree-sha1 = "960639a12ffd223ee463e93392aeb260fa325566" +uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5" +version = "0.5.8" + +[[WinRPM]] +deps = ["BinDeps", "Compat", "HTTPClient", "LibExpat", "Libdl", "Libz", "URIParser"] +git-tree-sha1 = "2a889d320f3b77d17c037f295859fe570133cfbf" +uuid = "c17dfb99-b4f7-5aad-8812-456da1ad7187" +version = "0.4.2" diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..9ebc032 --- /dev/null +++ b/Project.toml @@ -0,0 +1,24 @@ +name = "ParticleScattering" +uuid = "3532a73c-0c3e-5edc-89e3-d771336b0a2d" + +[deps] +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8" +LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/README.md b/README.md index 36dd8f2..88341c0 100644 --- a/README.md +++ b/README.md @@ -14,13 +14,16 @@ particle parameters for various design problems. ### Installation -ParticleScattering can be installed using `Pkg.add`. Currently, only Julia 0.6 is supported. +ParticleScattering for julia 0.7 can be installed using `Pkg.add`: ```julia Pkg.add("ParticleScattering") using ParticleScattering ``` +For julia 0.6, an older version of ParticleScattering can be installed manually +by cloning [release v0.0.4 from GitHub](https://github.com/bblankrot/ParticleScattering.jl/releases/tag/v0.0.4). + ### Community The easiest way to contribute is by opening issues! Of course, we'd be more than happy if you implement any fixes and send a PR. diff --git a/REQUIRE b/REQUIRE deleted file mode 100644 index db375bc..0000000 --- a/REQUIRE +++ /dev/null @@ -1,10 +0,0 @@ -julia 0.6 - -IterativeSolvers 0.5.0 -LinearMaps -Optim 0.14.0 -PyPlot -DataFrames -CSV -PGFPlotsX 0.2.1 -JLD diff --git a/appveyor.yml b/appveyor.yml index 08b9366..945608c 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,18 +1,11 @@ environment: matrix: - - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe" - PYTHONDIR: "use_conda" - PYTHON: "use_conda" - - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe" - PYTHONDIR: "use_conda" - PYTHON: "use_conda" + - julia_version: 0.7 + PYTHON: "Conda" -## uncomment the following lines to allow failures on nightly julia -## (tests will run but not make your overall status red) -#matrix: -# allow_failures: -# - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe" -# - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe" +platform: + - x86 # 32-bit + - x64 # 64-bit branches: only: @@ -26,24 +19,12 @@ notifications: on_build_status_changed: false install: - - ps: "[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12" -# If there's a newer build queued for the same PR, cancel this one - - ps: if ($env:APPVEYOR_PULL_REQUEST_NUMBER -and $env:APPVEYOR_BUILD_NUMBER -ne ((Invoke-RestMethod ` - https://ci.appveyor.com/api/projects/$env:APPVEYOR_ACCOUNT_NAME/$env:APPVEYOR_PROJECT_SLUG/history?recordsNumber=50).builds | ` - Where-Object pullRequestId -eq $env:APPVEYOR_PULL_REQUEST_NUMBER)[0].buildNumber) { ` - throw "There are newer queued builds for this pull request, failing early." } -# Download most recent Julia Windows binary - - ps: (new-object net.webclient).DownloadFile( - $env:JULIA_URL, - "C:\projects\julia-binary.exe") -# Run installer silently, output to C:\projects\julia - - C:\projects\julia-binary.exe /S /D=C:\projects\julia + - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1")) build_script: -# Need to convert from shallow to complete for Pkg.clone to work - - IF EXIST .git\shallow (git fetch --unshallow) - - C:\projects\julia\bin\julia -e "versioninfo(); - Pkg.clone(pwd(), \"ParticleScattering\"); Pkg.build(\"ParticleScattering\"); using Conda; Conda.add(\"basemap\");" + - echo "%JL_BUILD_SCRIPT%" + - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%" test_script: - - C:\projects\julia\bin\julia -e "Pkg.test(\"ParticleScattering\")" + - echo "%JL_TEST_SCRIPT%" + - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%" diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..cb41828 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,5 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "~0.22" diff --git a/docs/make.jl b/docs/make.jl index 2f210cd..2068fc3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,9 +1,11 @@ using Documenter, ParticleScattering -makedocs(format = :html, +makedocs( sitename = "ParticleScattering.jl", authors = "Boaz Blankrot", - linkcheck = !("skiplinks" in ARGS), + format = Documenter.HTML( + prettyurls = get(ENV, "CI", nothing) == "true" + ), pages = Any[ "Home" => "index.md", "Tutorials" => Any[ @@ -21,8 +23,4 @@ makedocs(format = :html, deploydocs( repo = "github.com/bblankrot/ParticleScattering.jl.git", - target = "build", - julia = "0.6", - deps = nothing, - make = nothing ) diff --git a/docs/src/assets/direct_vs_fmm0.png b/docs/src/assets/direct_vs_fmm0.png index e77f5ee..7af9e1b 100644 Binary files a/docs/src/assets/direct_vs_fmm0.png and b/docs/src/assets/direct_vs_fmm0.png differ diff --git a/docs/src/assets/fmm_tutorial_plot0.png b/docs/src/assets/fmm_tutorial_plot0.png index 17d07ab..5f77c87 100644 Binary files a/docs/src/assets/fmm_tutorial_plot0.png and b/docs/src/assets/fmm_tutorial_plot0.png differ diff --git a/docs/src/assets/optim_angle.png b/docs/src/assets/optim_angle.png index a3cca75..cc32112 100644 Binary files a/docs/src/assets/optim_angle.png and b/docs/src/assets/optim_angle.png differ diff --git a/docs/src/assets/optim_angle_conv.png b/docs/src/assets/optim_angle_conv.png index 49a38f0..2b5cd96 100644 Binary files a/docs/src/assets/optim_angle_conv.png and b/docs/src/assets/optim_angle_conv.png differ diff --git a/docs/src/assets/optim_radius_conv.png b/docs/src/assets/optim_radius_conv.png index e30aa04..bcd3d21 100644 Binary files a/docs/src/assets/optim_radius_conv.png and b/docs/src/assets/optim_radius_conv.png differ diff --git a/docs/src/incident_fields.md b/docs/src/incident_fields.md index ea17cb9..434b729 100644 --- a/docs/src/incident_fields.md +++ b/docs/src/incident_fields.md @@ -32,14 +32,14 @@ sources must lie completely outside all particles, and are contructed with `Curr ```julia using ParticleScattering, PyPlot λ = 1 -yc = linspace(-0.5λ, 0.5λ, 100) +yc = range(-0.5λ, stop=0.5λ, length=100) ui = CurrentSource(0, -0.5λ, 0, 0.5λ, cos.(π*yc/λ)) #now plot x_points = 100; y_points = 100 -x = linspace(-3λ, 3λ, x_points + 1) -y = linspace(-3λ, 3λ, y_points + 1) -xgrid = repmat(x', y_points + 1, 1) -ygrid = repmat(y, 1, x_points + 1) +x = range(-3λ, stop=3λ, length=x_points + 1) +y = range(-3λ, stop=3λ, length=y_points + 1) +xgrid = repeat(transpose(x), y_points + 1) +ygrid = repeat(y, 1, x_points + 1) points = cat(2, vec(xgrid[1:y_points, 1:x_points]) + 3λ/x_points, vec(ygrid[1:y_points, 1:x_points]) + 3λ/y_points) u = uinc(2π/λ, points, ui) diff --git a/docs/src/index.md b/docs/src/index.md index 6bd3ab7..5681705 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -7,23 +7,15 @@ particles. Provides the ability to optimize over particle parameters for various ## Installation -ParticleScattering can be installed using `Pkg.add`. Currently, only Julia 0.6 is supported. +ParticleScattering for julia 0.7 can be installed using `Pkg.add`: ```julia Pkg.add("ParticleScattering") using ParticleScattering ``` -which also installs the following dependencies: -```julia -IterativeSolvers -LinearMaps -Optim -PyPlot -DataFrames -CSV -PGFPlotsX -``` +For julia 0.6, an older version of ParticleScattering can be installed manually +by cloning [release v0.0.4 from GitHub](https://github.com/bblankrot/ParticleScattering.jl/releases/tag/v0.0.4). ## Getting Started diff --git a/docs/src/minimalNP.md b/docs/src/minimalNP.md index ad15ad0..584eec8 100644 --- a/docs/src/minimalNP.md +++ b/docs/src/minimalNP.md @@ -31,11 +31,7 @@ accuracy due to `N` is mostly irrelevant). Above a certain `N`, this error tends to decay as ``O(N^{-3})``, but with a multiplicative factor that is heavily dependent on the particle and wavelength. -With `minimumN`, we first guess a value and then use a binary search to find the -minimal `N` satisfying some error tolerance `tol`: - -```@docs -minimumN -``` +With [`minimumN`](@ref), we first guess a value and then use a binary search to find the +minimal `N` satisfying some error tolerance `tol`. ## Multipole Parameter P diff --git a/docs/src/tutorial1.md b/docs/src/tutorial1.md index 6d3f80b..ff14f65 100644 --- a/docs/src/tutorial1.md +++ b/docs/src/tutorial1.md @@ -50,7 +50,7 @@ appropriate method: ```julia #calculate field on the x-axis passing through the particle -points = [linspace(-0.5λ0, 0.5λ0, 200) zeros(200)] +points = [range(-0.5λ0, stop=0.5λ0, length=200) zeros(200)] u = calc_near_field(k0, kin, P, sp, points, pw) ``` diff --git a/docs/src/tutorial2.md b/docs/src/tutorial2.md index 8e361ff..2224c9e 100644 --- a/docs/src/tutorial2.md +++ b/docs/src/tutorial2.md @@ -66,11 +66,12 @@ converted to line sources, and are thus fully contained in the FMM grid. Calculating and plotting the near or far fields with FMM is just as in the [previous tutorial](@ref scattering_small_grid), except we must supply the -`FMMoptions` object: +`FMMoptions` object. We can also speed up the plotting by using the +`method = "recurrence"` option: ```julia plot_near_field(k0, kin, P, sp, PlaneWave(θ_i), opt = fmm_options, - border = [-12;12;-10;10], x_points = 480, y_points = 400) + border = [-12;12;-10;10], x_points = 480, y_points = 400, method = "recurrence") colorbar() ``` @@ -94,9 +95,9 @@ can arise in the direct solution of even the simplest scattering problems: k0 = 0.01 kin = 0.02 shapes = [squircle(1, 200)] -ids = [1;1] -centers = [0.0 0.0; 5.0 0.0] -phis = [0.0;0.0] +ids = ones(Int, 64) +centers = square_grid(8, 5) +phis = zeros(64) sp = ScatteringProblem(shapes, ids, centers, phis) Pmax = 15 ``` @@ -105,11 +106,11 @@ We solve this problem using the direct approach and with FMM, and then compare both the multipole coefficients ``\beta`` and the resulting potential densities: ```julia -betas = Array{Vector}(Pmax) -betas_FMM = Array{Vector}(Pmax) -inners = Array{Vector}(Pmax) -inners_FMM = Array{Vector}(Pmax) -fmmopts = ParticleScattering.FMMoptions(true, nx = 1, acc = 9) +betas = Array{Vector}(undef, Pmax) +betas_FMM = Array{Vector}(undef, Pmax) +inners = Array{Vector}(undef, Pmax) +inners_FMM = Array{Vector}(undef, Pmax) +fmmopts = ParticleScattering.FMMoptions(true, nx = 4, acc = 9) for P = 1:Pmax betas[P], inners[P] = solve_particle_scattering(k0, kin, P, sp, PlaneWave(); verbose = false) @@ -117,32 +118,33 @@ for P = 1:Pmax fmmopts; verbose = false) betas_FMM[P] = res[1] end - +import LinearAlgebra.norm errnorm(x,y) = norm(x-y)/norm(x) figure() subplot(2,1,1) -semilogy([errnorm(betas_FMM[i], betas[i]) for i = 1:Pmax]) +semilogy(1:Pmax, [errnorm(betas_FMM[i], betas[i]) for i = 1:Pmax]) +xlim([1;Pmax]); grid() ylabel("\$\\Delta \\beta\$") subplot(2,1,2) -semilogy([errnorm(inners_FMM[i], inners[i]) for i = 1:Pmax]) +semilogy(1:Pmax, [errnorm(inners_FMM[i], inners[i]) for i = 1:Pmax]) xlabel("\$ P \$") +xlim([1;Pmax]); grid() ylabel("\$ \\Delta \$" * " Potential Density") ``` -![direct_vs_fmm0](./assets/direct_vs_fmm0.png) + +```@raw html +

direct_vs_fmm0

+``` In both subplots, we see that increasing `P` actually leads to a decrease in accuracy (plotting the results separately also shows that the FMM results stay -the same, while the direct results blow up). This is due to two main reasons - +virtually constant, while the direct results blow up). This is due to two main reasons – conditioning of the system matrix, and the fact that high-order cylindrical harmonics are responsible for substantially greater potential densities than lower-order ones. Both of these are impacted by the number of particles as well -as the wavelength. +as the wavelength, but mitigated by the iterative solver used by the FMM solver. This ties in with [Choosing Minimal N and P](@ref minimalNP) – not only does increasing `P` far beyond that required for a certain error impact runtime, but can also increase the error in the solution. - -Of course, FMM was not really used here as `nx == 1` means both particles are in -the same FMM group, and the maintained accuracy is purely due to the iterative -system matrix solution used in `solve_particle_scattering_FMM`, GMRES. diff --git a/docs/src/tutorial_optim_angle.md b/docs/src/tutorial_optim_angle.md index f49c862..a18c9fb 100644 --- a/docs/src/tutorial_optim_angle.md +++ b/docs/src/tutorial_optim_angle.md @@ -1,4 +1,4 @@ -# Tutorial 3: Angle Optimization +# [Tutorial 3: Angle Optimization](@id tutorial3) In this tutorial, we build upon the [previous tutorial](@ref tutorial2) by optimizing the rotation angles of the particles (`φs`) to maximize the field @@ -40,19 +40,19 @@ where `φs0` is the starting point for the optimization method, and `points` are the locations at which we intend to maximize or minimize the field intensity. In this case, we want to optimize intensity at a small area around the origin. We now select the optimization method and select its options. -In most cases, this combination of BFGS with a backtracking line search will -yield accurate results in fast time; -other line searches that require re-evaluation of the gradient will be -significantly slower but may converge more accurately. The possible convergence +In most cases, BFGS with the default linesearch will yield accurate results in +fast time; optimizing with `adjoint = false` requires a backtracking linesearch +to reduce gradient evaluations. The possible convergence criteria are set by the bounds `f_tol`, `g_tol`, and `x_tol`, for a relative change in the function, gradient norm, or variables, respectively. In addition, we can set a maximum number of `iterations`. Verbosity of the output is set with `show_trace` and `extended_trace`. ```julia +import Optim optim_options = Optim.Options(f_tol = 1e-6, iterations = 100, store_trace = true, show_trace = true) -optim_method = Optim.BFGS(;linesearch = LineSearches.BackTracking()) +optim_method = Optim.BFGS() ``` We now run both minimization and maximization: @@ -71,41 +71,43 @@ separately with [`plot_near_field`](@ref) or compare them side by side with the following PyPlot code: ```julia -plts = Array{Any}(3) -plts[1] = plot_near_field(k0, kin, P, sp, pw, x_points = 100, y_points = 300, - opt = fmm_options, border = find_border(sp, points)) -plts[2] = plot_near_field(k0, kin, P, sp_min, pw, x_points = 100, y_points = 300, - opt = fmm_options, border = find_border(sp, points)) -plts[3] = plot_near_field(k0, kin, P, sp_max, pw, x_points = 100, y_points = 300, - opt = fmm_options, border = find_border(sp, points)) -close("all") - -fig, axs = subplots(ncols=3); msh = 0 +border = find_border(sp, points) +plts = Array{Any}(undef, 3) +for (i,spi) in enumerate([sp;sp_min;sp_max]) + global plts[i] = plot_near_field(k0, kin, P, spi, pw, x_points = 100, + y_points = 300, opt = fmm_options, border = border) +end + +using PyPlot +fig, axs = subplots(ncols=3, figsize=[4.5,4]) for (i, spi) in enumerate([sp;sp_min;sp_max]) - msh = axs[i][:pcolormesh](plts[i][2][1], plts[i][2][2], abs.(plts[i][2][3]), + global msh = axs[i].pcolormesh(plts[i][2][1], plts[i][2][2], abs.(plts[i][2][3]), vmin = 0, vmax = 3.4, cmap="viridis") draw_shapes(spi, ax = axs[i]) - axs[i][:set_aspect]("equal", adjustable = "box") - axs[i][:set_xlim]([border[1];border[2]]) - axs[i][:set_ylim]([border[3];border[4]]) - axs[i][:set_xticks]([-1,0,1]) - i > 1 && axs[i][:set_yticks]([]) + axs[i].set_aspect("equal", adjustable = "box") + axs[i].set_xlim([border[1];border[2]]) + axs[i].set_ylim([border[3];border[4]]) + axs[i].set_xticks([-1,0,1]) + i > 1 && axs[i].set_yticks([]) + axs[i].set_xlabel("x") end -subplots_adjust(left=0.05, right=0.8, top=0.98, bottom = 0.05, wspace = 0.1) -cbar_ax = fig[:add_axes]([0.85, 0.05, 0.05, 0.93]) -fig[:colorbar](msh, cax=cbar_ax) +axs[1].set_ylabel("y") +subplots_adjust(left=0.12, right=0.8, top=0.98, bottom = 0.12, wspace = 0.1) +cbar_ax = fig.add_axes([0.83, 0.12, 0.05, 0.86]) +fig.colorbar(msh, cax=cbar_ax) +cbar_ax.set_ylabel("|E\$_z\$ [V/m]") ``` ```@raw html
-optim_angle +optim_angle

``` From left to right, we see the electric field before optimization, after minimization, and after maximization. The field intensity at the origin is notably different in both optimization results, with minimization decreasing -the intensity by 95%, and maximization increasing it by over 700%. The convergence -of the optimization method for both examples can be plotted via: +the intensity at `points` by 77%, and maximization almost tripling it. +The convergence of the optimization method for both examples can be plotted via: ```julia iters = length(res_min.trace) @@ -115,20 +117,22 @@ iters2 = length(res_max.trace) fobj2 = -[res_max.trace[i].value for i=1:iters2] gobj2 = [res_max.trace[i].g_norm for i=1:iters2] -fig, axs = subplots(ncols=2, figsize=[7,5]) -axs[1][:semilogy](0:iters-1, fobj, linewidth=2) -axs[2][:semilogy](0:iters-1, gobj, linewidth=2) -axs[1][:semilogy](0:iters2-1, fobj2, linewidth=2, "--") -axs[2][:semilogy](0:iters2-1, gobj2, linewidth=2, "--") -axs[1][:legend](["\$f_\\mathrm{obj}\$ (min)"; +fig, axs = subplots(ncols=2, figsize=[6,3]) +axs[1].semilogy(0:iters-1, fobj, linewidth=2) +axs[2].semilogy(0:iters-1, gobj, linewidth=2) +axs[1].semilogy(0:iters2-1, fobj2, linewidth=2, "--") +axs[2].semilogy(0:iters2-1, gobj2, linewidth=2, "--") +axs[1].legend(["\$f_\\mathrm{obj}\$ (min)"; "\$f_\\mathrm{obj}\$ (max)"], loc="right") -axs[2][:legend](["\$\\|\\mathbf{g}_\\mathrm{obj}\\|\$ (min)"; - "\$\\|\\mathbf{g}_\\mathrm{obj}\\|\$ (max)"], loc="best") -axs[1][:set_xlabel]("Iteration") -axs[2][:set_xlabel]("Iteration") -axs[1][:set_ylim](ymax=40) +axs[2].legend(["\$\\Vert \\nabla f_\\mathrm{obj}\\Vert\$ (min)"; + "\$\\Vert \\nabla f_\\mathrm{obj}\\Vert\$ (max)"], loc="best") +axs[1].set_xlabel("Iteration") +axs[2].set_xlabel("Iteration") +axs[1].set_ylim(ymax=40) +axs[1].grid(); axs[2].grid() +subplots_adjust(top=0.99, right=0.99, left=0.07, bottom=0.15) ``` ```@raw html -

optim_angle_conv

+

optim_angle_conv

``` diff --git a/docs/src/tutorial_optim_angle_pwr.md b/docs/src/tutorial_optim_angle_pwr.md new file mode 100644 index 0000000..a327854 --- /dev/null +++ b/docs/src/tutorial_optim_angle_pwr.md @@ -0,0 +1,7 @@ +# Tutorial 5: Angle Optimization for Power Flow + +In the [first angle optimization tutorial](@ref tutorial3), we optimized the rotation angles `φs` such that the field intensity was minimized or maximized around the center of the structure. In some cases, however, we wish to optimize electromagnetic power flow ``P`` through an arc ``\ell``: + +```math +P = \frac{1}{2} \Re \int (\mathbf{E} \times \mathbf{H}^*) \cdot \hat{\mathbf{n}} \, \mathrm{d} \ell +``` diff --git a/docs/src/tutorial_optim_radius.md b/docs/src/tutorial_optim_radius.md index d9caab5..b1165dd 100644 --- a/docs/src/tutorial_optim_radius.md +++ b/docs/src/tutorial_optim_radius.md @@ -26,8 +26,7 @@ fmm_options = FMMoptions(true, acc = 6, dx = 2a) [`optimize_radius`](@ref) not only allows us to optimize all of the radii simultaneously, but also to assign several particles the same `id`, which can be useful when the target radii are expected to have symmetry of some type. -Here we shall assume symmetry with respect to the ``x``-axis (horizontal line -of symmetry) with `uniqueind`: +Here we shall assume symmetry with respect to the ``x``-axis with `uniqueind`: ```julia # let's impose symmetry wrt x-axis @@ -55,8 +54,7 @@ points = [4a 0.0] r_max = (0.4*a)*ones(J) r_min = (1e-3*a)*ones(J) rs0 = (0.25*a)*ones(J) -assert(verify_min_distance([CircleParams(r_max[i]) for i = 1:J], - centers, ids, points)) +@assert verify_min_distance(CircleParams.(r_max), centers, ids, points) ``` The optimization process is initiated by running: @@ -71,14 +69,14 @@ With the optimization process complete, we can plot the electric field with the initial and optimized radii: ```julia -sp1 = ScatteringProblem([CircleParams(rs0[i]) for i = 1:J], ids, centers, φs) +sp1 = ScatteringProblem(CircleParams.(rs0), ids, centers, φs) plot_near_field(k0, kin, P, sp1, pw, x_points = 150, y_points = 150, opt = fmm_options, border = 0.9*[-1;1;-1;1], normalize = a) colorbar() clim([0;2.5]) xlabel("x/a") ylabel("y/a") -sp2 = ScatteringProblem([CircleParams(rs[i]) for i = 1:J], ids, centers, φs) +sp2 = ScatteringProblem(CircleParams.(rs), ids, centers, φs) plot_near_field(k0, kin, P, sp2, pw, x_points = 150, y_points = 150, opt = fmm_options, border = 0.9*[-1;1;-1;1], normalize = a) colorbar() diff --git a/examples/sim0_min_NP.jl b/examples/sim0_min_NP.jl index c6a8ad5..d7d3299 100644 --- a/examples/sim0_min_NP.jl +++ b/examples/sim0_min_NP.jl @@ -1,7 +1,12 @@ ###################################### -## comparing analytical (fictitious sources) and computed fields -using PyPlot,ParticleScattering +## comparing analytical (fictitious sources) and computed fields. minimum N +## routine runs on multiple cores and uses a lot of memory -- reduce max Nvec5, +## Nvecs if necessary +using PyPlot,Distributed,SpecialFunctions import JLD +addprocs(4) +@everywhere using ParticleScattering +output_dir = homedir() function unique_subset(N,v,rev = true) all(sort(v, rev = rev) .== v) || error("Only handles sorted vectors.") @@ -18,67 +23,10 @@ function unique_subset(N,v,rev = true) end N_u,v_u end - -k0 = 10.0 -kin = 1.5*k0 -l0 = 2π/k0 -a1 = 0.3l0 -a2 = 0.1l0 -N_points = 20_000 -θ_i = 0.0 - -myshapefun5(N) = rounded_star(a1,a2,5,N) -myshapefun_squircle(N) = squircle(a1+0.5a2,N) - -shape_functions = [myshapefun5; myshapefun_squircle] - -if !isfile(dirname(@__FILE__) * "/mindata.jld") - Nvec5 = [collect(20:10:90);collect(100:20:6000)] - Nvecs = [collect(20:10:90);collect(100:20:6000)] - errNvec5 = Array{Float64}(length(Nvec5)) - errNvecs = Array{Float64}(length(Nvecs)) - - s = myshapefun5(400) #just for radius - err_points = [s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)] - E_ana = (0.25im*besselh(0,1,k0*s.R))*ones(Complex{Float64},N_points) - E_comp = Array{Complex{Float64}}(length(E_ana)) - - for i in eachindex(Nvec5) - errNvec5[i] = ParticleScattering.minimumN_helper(Nvec5[i], k0, kin, - myshapefun5, err_points, E_comp, E_ana) - display("Nvec5: done with $i/$(length(Nvec5))") - end - - s = myshapefun_squircle(400) #just for radius - err_points = [s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)] - E_ana = (0.25im*besselh(0,1,k0*s.R))*ones(Complex{Float64},N_points) - E_comp = Array{Complex{Float64}}(length(E_ana)) - - for i in eachindex(Nvecs) - errNvecs[i] = ParticleScattering.minimumN_helper(Nvecs[i], k0, kin, - myshapefun_squircle, err_points, E_comp, E_ana) - display("Nvecs: done with $i/$(length(Nvecs))") - end - N5, errN5 = unique_subset(Nvec5, errNvec5) - Ns, errNs = unique_subset(Nvecs, errNvecs) - JLD.@save (dirname(@__FILE__) * "/mindata.jld") k0 kin l0 a1 a2 N5 errN5 Ns errNs -else - JLD.@load (dirname(@__FILE__) * "/mindata.jld") -end - - -###############################3 - -inds = findfirst(errNs .< 5e-10) -ind5 = findfirst(errN5 .< 5e-10) -N5 = N5[1:ind5] -Ns = Ns[1:inds] -errN5 = errN5[1:ind5] -errNs = errNs[1:inds] - -# here binary search is suboptimal since we know that new P is "close" to old one, frequently P or P+1 +# binary search here is suboptimal since we know that new P is "close" to old one, +# frequently P or P+1 function findMinP(N, errN, shapefun, N_points, k0, kin; P_last = 1, P_max = 100) - E_multipole = Array{Complex{Float64}}(N_points) + E_multipole = Array{Complex{Float64}}(undef, N_points) errP = zeros(Float64, length(errN)) Pmin = zeros(Int, length(errN)) for iN in eachindex(errN) @@ -97,86 +45,91 @@ function findMinP(N, errN, shapefun, N_points, k0, kin; P_last = 1, P_max = 100) P_last = P break elseif P == P_max - warn("Failed to find P for iN = $iN") + @warn("Failed to find P for iN = $iN") return Pmin, errP end end - display("findMinP: done with iN=$iN/$(length(errN)), matched P=$(Pmin[iN]) with $(errP[iN]) ≦ $(errN[iN])") + display("findMinP: done with iN=$iN/$(length(errN)), matched P=$(Pmin[iN]) with $(errP[iN]) ≤ $(errN[iN])") end Pmin, errP end -Ps,errPs = findMinP(Ns, errNs, myshapefun_squircle, N_points, k0, kin) -JLD.@save dirname(@__FILE__) * "/mindata_s_new.jld" Ns errNs Ps errPs - -P5,errP5 = findMinP(N5, errN5, myshapefun5, N_points, k0, kin) -JLD.@save dirname(@__FILE__) * "/mindata_5_new.jld" N5 errN5 P5 errP5 - +@everywhere begin + k0 = 10.0 + kin = 1.5*k0 + l0 = 2π/k0 + a1 = 0.3l0 + a2 = 0.1l0 + N_points = 20_000 + myshapefun5(N) = rounded_star(a1,a2,5,N) + myshapefun_squircle(N) = squircle(a1+0.5a2,N) + Nvec5 = unique(round.(Int, 10 .^ range(log10(20), stop=log10(5000), length=200))) + Nvecs = unique(round.(Int, 10 .^ range(log10(20), stop=log10(5000), length=200))) +end -##################################################3 -#plot with pgfplots -import PGFPlotsX; const pgf = PGFPlotsX +dt_Nvec5 = @elapsed begin + s = myshapefun5(400) #just for radius + err_points = [s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)] + E_ana = (0.25im*besselh(0,1,k0*s.R))*ones(Complex{Float64},N_points) + E_comp = Array{ComplexF64}(undef, length(E_ana)) -JLD.@load dirname(@__FILE__) * "/mindata_5_new.jld" -JLD.@load dirname(@__FILE__) * "/mindata_s_new.jld" + innerfunc1 = function(i) + ParticleScattering.minimumN_helper(Nvec5[i], k0, kin, + myshapefun5, err_points, E_comp, E_ana) + end + errNvec5 = pmap(innerfunc1, eachindex(Nvec5), + on_error = ex->(isa(ex, OutOfMemoryError) ? NaN : rethrow(ex))) + #or + #errNvec5 = pmap(innerfunc1, eachindex(Nvec5)) +end +display("finished calculating minimum N for rounded star in $dt_Nvec5 seconds") +any(isnan.(errNvec5)) && @warn "encountered OutOfMemoryError, some values are NaN" -pgf.@pgf begin - N_plot5 = pgf.Plot({blue, dashdotdotted, no_markers, thick}, - pgf.Coordinates(errN5,N5)) - N_plot = pgf.Plot({red, dotted, no_markers, thick}, - pgf.Coordinates(errNs,Ns)) - P_plot5 = pgf.Plot({"green!50!black", no_markers, thick}, - pgf.Coordinates(errN5,P5)) - P_plot = pgf.Plot({black, dashed, no_markers, thick}, - pgf.Coordinates(errNs,Ps)) - ax = pgf.Axis( - { - xmin = 5e-10, - xmax = 1e-1, - xlabel = "\$\\Delta u\$", - xmode = "log", - ymode = "log", - grid = "both", - legend_pos = "north east", - legend_style = "font = \\footnotesize", - legend_cell_align = "left" - }, - [N_plot5, N_plot, P_plot5, P_plot], - pgf.Legend(["\$N_{\\mathrm{min}} \\, \\mathrm{(star)}\$", - "\$N_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$", - "\$P_{\\mathrm{min}} \\, \\mathrm{(star)}\$", - "\$P_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$"])) +dt_Nvecs = @elapsed begin + s = myshapefun_squircle(400) #just for radius + err_points = [s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)] + E_ana = (0.25im*besselh(0,1,k0*s.R))*ones(Complex{Float64},N_points) + E_comp = Array{ComplexF64}(undef, length(E_ana)) + innerfunc2 = function(i) + ParticleScattering.minimumN_helper(Nvecs[i], k0, kin, + myshapefun_squircle, err_points, E_comp, E_ana) + end + errNvecs = pmap(innerfunc2, eachindex(Nvecs), + on_error = ex->(isa(ex, OutOfMemoryError) ? NaN : rethrow(ex))) + #or + #errNvecs = pmap(innerfunc2, eachindex(Nvecs)) end -pgf.save(dirname(@__FILE__) * "/tikz/minNP.tex", ax ,include_preamble = false) +display("finished calculating minimum N for squircle in $dt_Nvecs seconds") +any(isnan.(errNvecs)) && @warn "encountered OutOfMemoryError, some values are NaN" -shape = myshapefun5(200) -pgf_shape = pgf.Coordinates([shape.ft[:,1];shape.ft[1,1]], - [shape.ft[:,2];shape.ft[1,2]]) -pgf.@pgf begin - ax2 = pgf.Axis({axis_equal, ticks = "none"}, pgf.Plot({thick, black, - no_markers, fill = "black!20"}, pgf_shape)) - Rd = shape.R - push!(ax2, "\\addplot [black, dashed, thick, domain=0:2*pi,samples=100]({$(shape.R)*cos(deg(x))},{$(shape.R)*sin(deg(x))});") - push!(ax2, "\\node at (axis cs: -0.28,-0.15) {\$D\$};") - push!(ax2, "\\node at (axis cs: 0.25,-0.25) {\$k_0\$};") - push!(ax2, "\\node at (axis cs: 0.0,0.0) {\$k_1\$};") - push!(ax2, "\\node at (axis cs: 0.15,0.1) {\$\\partial \\Omega_1\$};") -end -pgf.save(dirname(@__FILE__) * "/tikz/inclusion.tex", ax2 ,include_preamble = false) +N5, errN5 = unique_subset(Nvec5, errNvec5) +Ns, errNs = unique_subset(Nvecs, errNvecs) +inds = something(findfirst(errNs .< 5e-10), length(errNs)) +ind5 = something(findfirst(errN5 .< 5e-10), length(errN5)) +N5 = N5[1:ind5] +Ns = Ns[1:inds] +errN5 = errN5[1:ind5] +errNs = errNs[1:inds] +P5,errP5 = findMinP(N5, errN5, myshapefun5, N_points, k0, kin) +Ps,errPs = findMinP(Ns, errNs, myshapefun_squircle, N_points, k0, kin) -shape = myshapefun_squircle(200) -pgf_shape = pgf.Coordinates([shape.ft[:,1];shape.ft[1,1]], - [shape.ft[:,2];shape.ft[1,2]]) -pgf.@pgf begin - ax3 = pgf.Axis({axis_equal, ticks = "none"}, pgf.Plot({thick, black, - no_markers, fill = "black!20", label="bla"}, pgf_shape)) - Rd = shape.R - push!(ax3, "\\addplot [black, dashed, thick, domain=0:2*pi,samples=100]({$(shape.R)*cos(deg(x))},{$(shape.R)*sin(deg(x))});") - push!(ax3, "\\node at (axis cs: -0.28,-0.15) {\$D\$};") - push!(ax3, "\\node at (axis cs: 0.25,-0.25) {\$k_0\$};") - push!(ax3, "\\node at (axis cs: 0.0,0.0) {\$k_1\$};") - push!(ax3, "\\node at (axis cs: 0.15,0.1) {\$\\partial \\Omega_1\$};") -end -pgf.save(dirname(@__FILE__) * "/tikz/inclusion2.tex", ax3 ,include_preamble = false) +JLD.@save joinpath(output_dir, "mindata.jld") k0 kin l0 a1 a2 N5 errN5 Ns errNs Ps errPs P5 errP5 + +################################################## +#plot +fig = figure(figsize=[3.5, 2.5]) +loglog(errN5, N5, "b-.", linewidth=1, + label="\$N_{\\mathrm{min}} \\, \\mathrm{(star)}\$") +loglog(errNs, Ns, "r:", linewidth=1, + label="\$N_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$") +loglog(errN5, P5, "tab:green", linewidth=1, + label="\$P_{\\mathrm{min}} \\, \\mathrm{(star)}\$") +loglog(errNs, Ps, "k--", linewidth=1, + label="\$P_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$") +xlim([5e-10;1e-1]) +xlabel("\$\\Delta \\, u\$") +tick_params(which="both", direction="in") +legend(fontsize=8) +grid() diff --git a/examples/sim1_complexity.jl b/examples/sim1_complexity.jl index 9c17a83..bc6a512 100644 --- a/examples/sim1_complexity.jl +++ b/examples/sim1_complexity.jl @@ -1,14 +1,17 @@ # Here is a simulation of solution time as a function of M - number of shapes (or sqrt of number of shapes) using ParticleScattering, IterativeSolvers, PyPlot - +import JLD +import LinearMaps: LinearMap +import Statistics: mean +output_dir = homedir() #loop definitions sqrtM_vec = collect(5:30); M_vec = sqrtM_vec.^2 trials = 3 simlen = length(M_vec) -res_vec = Array{Float64}(simlen,trials) -iter_vec = Array{Float64}(simlen,trials) -mvp_vec = Array{Float64}(simlen,trials) -setup_vec = Array{Float64}(simlen,trials) +res_vec = Array{Float64}(undef, simlen, trials) +iter_vec = Array{Float64}(undef, simlen, trials) +mvp_vec = Array{Float64}(undef, simlen, trials) +setup_vec = Array{Float64}(undef, simlen, trials) #variables k0 = 10.0 @@ -27,77 +30,83 @@ P = 10; errP = 9.57e-7 # P,errP = minimumP(k0, kin, shapes[1], tol = tol, N_points = 20_000, # P_min = 1, P_max = 120) -tic() scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, [1]) α = Complex{Float64}[exp(1.0im*p*(pi/2-θ_i)) for p=-P:P] -dt0 = toq() -for i=1:simlen-1 - tic() - scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, [1]) - α = Complex{Float64}[exp(1.0im*p*(pi/2-θ_i)) for p=-P:P] - dt0 += toq() +dt0 = @elapsed for i=1:trials + global scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, [1]) + global α = Complex{Float64}[exp(1.0im*p*(pi/2-θ_i)) for p=-P:P] end -dt0 /= simlen +dt0 /= trials -for is = 1:simlen, it = 1:trials +function time_FMM(is, it) #compute shape variables - begin #setup - sqrtM = sqrtM_vec[is] - M = sqrtM^2 - centers = square_grid(sqrtM, dist) - φs = rand(M) - ids = ones(Int, M) - opt = FMMoptions(true, acc = Int(-log10(tol)), nx = div(sqrtM,2), method="pre") - end - tic() - (groups, boxSize) = divideSpace(centers, opt) - (P2, Q) = FMMtruncation(opt.acc, boxSize, k0) - mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri=true) + sqrtM = sqrtM_vec[is] + M = sqrtM^2 + centers = square_grid(sqrtM, dist) + φs = rand(M) + ids = ones(Int, M) + opt = FMMoptions(true, acc = Int(-log10(tol)), nx = div(sqrtM,2), method="pre") + + setup_vec[is,it] = @elapsed begin + (groups, boxSize) = divideSpace(centers, opt) + (P2, Q) = FMMtruncation(opt.acc, boxSize, k0) + mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri=true) - #construct rhs - rhs = repeat(α,outer=[M]) - for ic = 1:M - rng = (ic-1)*(2*P+1) + (1:2*P+1) - if φs[ic] == 0.0 - rhs[rng] = scatteringMatrices[ids[ic]]*α - else - #rotate without matrix - ParticleScattering.rotateMultipole!(view(rhs,rng),-φs[ic],P) - rhs[rng] = scatteringMatrices[ids[ic]]*rhs[rng] - ParticleScattering.rotateMultipole!(view(rhs,rng),φs[ic],P) + #construct rhs + rhs = repeat(α, M) + for ic = 1:M + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) + if φs[ic] == 0.0 + rhs[rng] = scatteringMatrices[ids[ic]]*α + else + #rotate without matrix + ParticleScattering.rotateMultipole!(view(rhs,rng),-φs[ic],P) + rhs[rng] = scatteringMatrices[ids[ic]]*rhs[rng] + ParticleScattering.rotateMultipole!(view(rhs,rng),φs[ic],P) + end + #phase shift added to move cylinder coords + phase = exp(1.0im*k0*(cos(θ_i)*centers[ic,1] + sin(θ_i)*centers[ic,2])) + rhs[rng] .*= phase end - #phase shift added to move cylinder coords - phase = exp(1.0im*k0*(cos(θ_i)*centers[ic,1] + sin(θ_i)*centers[ic,2])) - rhs[rng] .*= phase - end - pre_agg_buffer = zeros(Complex{Float64},Q,length(groups)) - trans_buffer = Array{Complex{Float64}}(Q) + pre_agg_buffer = zeros(Complex{Float64},Q,length(groups)) + trans_buffer = Array{Complex{Float64}}(undef, Q) - MVP = LinearMaps.LinearMap{eltype(rhs)}((output_, x_) -> - ParticleScattering.FMM_mainMVP_pre!(output_, x_, scatteringMatrices, - φs, ids, P, mFMM, pre_agg_buffer, trans_buffer), - M*(2*P+1), M*(2*P+1), ismutating = true) - x = zero(rhs) - setup_vec[is,it] = toq() + MVP = LinearMap{eltype(rhs)}((output_, x_) -> + ParticleScattering.FMM_mainMVP_pre!(output_, x_, scatteringMatrices, + φs, ids, P, mFMM, pre_agg_buffer, trans_buffer), + M*(2*P+1), M*(2*P+1), ismutating = true) + x = zero(rhs) + end - tic() - x,ch = gmres!(x, MVP, rhs, restart = M*(2*P+1), tol = opt.tol, log = true, - initially_zero = true) #no restart, preconditioning - res_vec[is,it] = toq() - tic() - rhs[:] = MVP*x - mvp_vec[is,it] = toq() + res_vec[is,it] = @elapsed begin + x,ch = gmres!(x, MVP, rhs, restart = M*(2*P+1), tol = opt.tol, + log = true, initially_zero = true) + end + mvp_vec[is,it] = @elapsed begin + rhs[:] = MVP*x + end iter_vec[is,it] = ch.iters end - -using JLD; @save "sim1_complexity.jld" res_vec, mvp_vec +#warmup +for is = 1:5:simlen + time_FMM(is, 1) +end +display("starting main benchmark...") +#benchmark +for is = 1:simlen, it = 1:trials + time_FMM(is, it) +end #average over all simulations -res_vec = vec(mean(res_vec,2)) -iter_vec = vec(mean(iter_vec,2)) -mvp_vec = vec(mean(mvp_vec,2)) -setup_vec = vec(mean(setup_vec,2)) +res_vec = vec(mean(res_vec, dims=2)) +iter_vec = vec(mean(iter_vec, dims=2)) +mvp_vec = vec(mean(mvp_vec, dims=2)) +setup_vec = vec(mean(setup_vec, dims=2)) +JLD.@save(joinpath(output_dir, "complexity.jld"), M_vec, res_vec, mvp_vec, setup_vec) + +######################## +linreg(x, y) = hcat(fill!(similar(x), 1), x) \ y a_total,b_total = linreg(log10.(M_vec), log10.(res_vec)) res_ana = (10^a_total)*(M_vec.^b_total) a_mvp,b_mvp = linreg(log10.(M_vec), log10.(mvp_vec)) @@ -114,35 +123,3 @@ legend(("Elapsed time (Sol.)", @sprintf("\$%fM^{%.2f}\$", 10^a_total, b_total), "Elapsed time (MVP.)", @sprintf("\$%fM^{%.2f}\$", 10^a_mvp, b_mvp), "Elapsed time (Setup)", @sprintf("\$%fM^{%.2f}\$", 10^a_setup, b_setup)), loc = "best") xlabel("Number of Scatterers") - -### plot with pgfplots -import PGFPlotsX; const pgf = PGFPlotsX -pgf.@pgf begin - ax = pgf.Axis({xlabel = "Number of scatterers", - ylabel = "\$ \\mathrm{Run time} \\ [\\mathrm{s}]\$", - xmode = "linear", - ymode = "log", - width = "10cm", - legend_pos = "north west", - legend_style = "font = \\footnotesize", - legend_cell_align = "left"}) - push!(ax, pgf.Plot({blue, "only marks", mark = "*"}, - pgf.Coordinates(M_vec, res_vec))) - tmp_total = floor(log10(10^a_total)) - push!(ax, pgf.Plot({blue, thick, no_markers}, - pgf.Coordinates(M_vec, res_ana))) - push!(ax, pgf.Plot({red, only_marks, mark = "triangle*", - mark_options = {fill = "red"}}, - pgf.Coordinates(M_vec, mvp_vec))) - tmp_mvp = floor(log10(10^a_mvp)) - push!(ax, pgf.Plot({red, thick, dashed, no_markers}, - pgf.Coordinates(M_vec, mvp_ana))) - push!(ax, pgf.Legend([ - "Elapsed time (solution)", - @sprintf("\$%.1f \\cdot 10^{%d} \\cdot M^{%.2f}\$", - 10^a_total/10^tmp_total, tmp_total, b_total), - "Elapsed time (MVP)", - @sprintf("\$%.1f \\cdot 10^{%d} \\cdot M^{%.2f}\$", - 10^a_mvp/10^tmp_mvp, tmp_mvp, b_mvp)])) -end -pgf.save(dirname(@__FILE__) * "/sim1.tex", ax ,include_preamble = false) diff --git a/examples/sim2_LuneburgOptim.jl b/examples/sim2_LuneburgOptim.jl index 9142853..612edba 100644 --- a/examples/sim2_LuneburgOptim.jl +++ b/examples/sim2_LuneburgOptim.jl @@ -12,7 +12,7 @@ kin = k0*sqrt(er) N_cells = Int(round(2*R_lens/a_lens)) centers, ids_lnbrg, rs_lnbrg = luneburg_grid(R_lens, N_cells, er) φs = zeros(Float64, length(ids_lnbrg)) -θ_i = 0.0 +pw = PlaneWave(0.0) P = 5 fmm_options = FMMoptions(true, acc = 6, dx = 2*a_lens, method = "pre") @@ -28,40 +28,31 @@ r_min = (a_lens*1e-3)*ones(size(centers,1)) rs0 = (0.25*a_lens)*ones(size(centers,1)) ids_max = collect(1:length(rs0)) -test_max = optimize_radius(rs0, r_min, r_max, points, ids_max, P, PlaneWave(θ_i), k0, kin, #precompile +test_max = optimize_radius(rs0, r_min, r_max, points, ids_max, P, pw, k0, kin, #precompile centers, fmm_options, optim_options, minimize = false) -tic() -test_max = optimize_radius(rs0, r_min, r_max, points, ids_max, P, PlaneWave(θ_i), k0, kin, +optim_time = @elapsed begin + test_max = optimize_radius(rs0, r_min, r_max, points, ids_max, P, pw, k0, kin, centers, fmm_options, optim_options, minimize = false) -optim_time = toq() +end rs_max = test_max.minimizer # plot near fields -filename1 = output_dir * "/opt_r_luneburg.tex" -filename2 = output_dir * "/opt_r_max.tex" -filename3 = output_dir * "/opt_r_0.tex" border = (R_lens + a_lens)*[-1;1;-1;1] sp1 = ScatteringProblem([CircleParams(rs_lnbrg[i]) for i in eachindex(rs_lnbrg)], ids_lnbrg, centers, φs) -Ez1 = plot_near_field(k0, kin, P, sp1, PlaneWave(θ_i), x_points = 150, y_points = 150, +Ez1 = plot_near_field(k0, kin, P, sp1, pw, x_points = 300, y_points = 300, opt = fmm_options, border = border) -plot_near_field_pgf(filename1, k0, kin, P, sp1, PlaneWave(θ_i); opt = fmm_options, - x_points = 201, y_points = 201, border = border) sp2 = ScatteringProblem([CircleParams(rs_max[i]) for i in eachindex(rs_max)], ids_max, centers, φs) -Ez2 = plot_near_field(k0, kin, P, sp2, PlaneWave(θ_i), x_points = 150, y_points = 150, +Ez2 = plot_near_field(k0, kin, P, sp2, pw, x_points = 300, y_points = 300, opt = fmm_options, border = border) -plot_near_field_pgf(filename2, k0, kin, P, sp2, PlaneWave(θ_i); opt = fmm_options, - x_points = 201, y_points = 201, border = border) sp3 = ScatteringProblem([CircleParams(rs0[i]) for i in eachindex(rs0)], collect(1:length(rs0)), centers, φs) -Ez3 = plot_near_field(k0, kin, P, sp3, PlaneWave(θ_i), x_points = 150, y_points = 150, +Ez3 = plot_near_field(k0, kin, P, sp3, pw, x_points = 300, y_points = 300, opt = fmm_options, border = border) -plot_near_field_pgf(filename3, k0, kin, P, sp3, PlaneWave(θ_i); opt = fmm_options, - x_points = 201, y_points = 201, border = border) #plot convergence inner_iters = length(test_max.trace) @@ -70,43 +61,20 @@ fobj = -[test_max.trace[i].value for i=1:inner_iters] gobj = [test_max.trace[i].g_norm for i=1:inner_iters] rng = iters .== 0 -test_max_trace = test_max.trace trace_of_r = [test_max.trace[i].metadata["x"] for i=1:inner_iters] -JLD.@save output_dir * "/luneburg_optim.jld" - -# figure() -# plot(0:inner_iters-1, fobj) -# plot(0:inner_iters-1, gobj) -# plot((0:inner_iters-1)[rng], fobj[rng],"*") -# plot((0:inner_iters-1)[rng], gobj[rng],"*") - -import PGFPlotsX; const pgf = PGFPlotsX -pgf.@pgf begin - fobj_plot = pgf.Plot({blue, thick, no_markers}, - pgf.Coordinates(0:inner_iters-1, fobj)) - gobj_plot = pgf.Plot({red, thick, dashed, no_markers}, - pgf.Coordinates(0:inner_iters-1, gobj)) - fobj_outer = pgf.Plot({blue, only_marks, mark = "*", - mark_options = {fill = "blue"}}, - pgf.Coordinates((0:inner_iters-1)[rng], fobj[rng])) - gobj_outer = pgf.Plot({red, only_marks, mark = "triangle*", - mark_options = {fill = "red"}}, - pgf.Coordinates((0:inner_iters-1)[rng], gobj[rng])) - ax = pgf.Axis({ xmin = 0, - width = "6cm", - xlabel = "Iterations", - legend_pos = "north east", - legend_style = "font = \\footnotesize", - legend_cell_align = "left"}, - fobj_plot, gobj_plot, fobj_outer, gobj_outer) - push!(ax, pgf.Legend(["\$f_{\\mathrm{obj}}\$"; - "\$\\|\\mathbf{g}_{\\mathrm{obj}}\\|_{\\infty}\$"])) - ax -end -pgf.save(output_dir * "/opt_r_conv.tex", ax ,include_preamble = false) +JLD.@save(joinpath(output_dir, "luneburg_optim.jld"), sp1, sp2, sp3,Ez1, Ez2, Ez3, + inner_iters, iters, fobj, gobj, rng, optim_time, rs_max, border) + +figure() +plot(0:inner_iters-1, fobj, "b", label="\$f_{\\mathrm{obj}}\$") +plot(0:inner_iters-1, log10.(gobj), "r--", label="\$\\Vert \\mathbf{g}_{\\mathrm{obj}}\\Vert_{\\infty}\$") +plot((0:inner_iters-1)[rng], fobj[rng],"bo") +plot((0:inner_iters-1)[rng], gobj[rng],"r^") +xlabel("Iterations") +legend(loc="best") ################ Testing with symmetry ###################### -assert(length(ids_max)==size(centers,1)) +@assert length(ids_max)==size(centers,1) centers_abs = centers[:,1] + 1im*abs.(centers[:,2]) ids_sym, centers_abs = ParticleScattering.uniqueind(centers_abs) J = length(centers_abs) @@ -114,27 +82,20 @@ r_max = (a_lens/1.15/2)*ones(J) r_min = (a_lens*1e-3)*ones(J) rs0 = (0.25*a_lens)*ones(J) -tic() -test_max_sym = optimize_radius(rs0, r_min, r_max, points, ids_sym, P, PlaneWave(θ_i), k0, kin, +sym_time = @elapsed begin + test_max_sym = optimize_radius(rs0, r_min, r_max, points, ids_sym, P, pw, k0, kin, centers, fmm_options, optim_options, minimize = false) -sym_time = toq() -rs_sym = test_max_sym.minimizer -JLD.@save output_dir * "/luneburg_optim_sym.jld" test_max_sym sym_time +end +JLD.@save joinpath(output_dir, "luneburg_optim_sym.jld") test_max_sym sym_time +rs_sym = test_max_sym.minimizer sp4 = ScatteringProblem([CircleParams(rs_sym[i]) for i in eachindex(rs_sym)], ids_sym, centers, φs) -Ez4 = plot_near_field(k0, kin, P, sp4, PlaneWave(θ_i), x_points = 150, y_points = 150, - opt = fmm_options, border = border) -pw = PlaneWave(θ_i) +Ez4 = plot_near_field(k0, kin, P, sp4, pw, x_points = 300, y_points = 300, + opt = fmm_options, border = border, method = "recurrence") + u1 = calc_near_field(k0, kin, 7, sp1, points, pw; opt = fmm_options) u2 = calc_near_field(k0, kin, 7, sp2, points, pw; opt = fmm_options) u3 = calc_near_field(k0, kin, 7, sp3, points, pw; opt = fmm_options) u4 = calc_near_field(k0, kin, 7, sp4, points, pw; opt = fmm_options) abs.([u1[1];u2[1];u3[1];u4[1]]) - -##################################### -# selfconsistent err P calculation -Ez_4 = calc_near_field(k0, kin, 4, sp1, points, pw; opt = fmm_options) -Ez_5 = calc_near_field(k0, kin, 5, sp1, points, pw; opt = fmm_options) -Ez_6 = calc_near_field(k0, kin, 6, sp1, points, pw; opt = fmm_options) -Ez_7 = calc_near_field(k0, kin, 7, sp1, points, pw; opt = fmm_options) diff --git a/examples/sim3_AngleOpt.jl b/examples/sim3_AngleOpt.jl index 6f428c8..b013ac5 100644 --- a/examples/sim3_AngleOpt.jl +++ b/examples/sim3_AngleOpt.jl @@ -1,6 +1,6 @@ using PyPlot, ParticleScattering -import JLD, Optim, PGFPlotsX; const pgf = PGFPlotsX - +import JLD, Optim, LineSearches +input_dir = homedir() output_dir = homedir() Ns = 100 k0 = 2π @@ -8,18 +8,17 @@ kin = 3*k0 l0 = 2π/k0 a1=0.3*l0; a2=0.1*l0; a3=5; dmin = R_multipole*2*(a1+a2) -θ_i = 0.5π +θ_i = 0.5π; ui = PlaneWave(θ_i) -size_factor = 7 -width = size_factor*3*l0 -height = size_factor*l0 +width = 21l0 +height = 7l0 myshapefun(N) = rounded_star(a1,a2,a3,N) -points = [linspace(0.0,width,10) height*ones(10)] +points = [range(0.0, stop=width, length=20) height*ones(20)] -if !isfile(dirname(@__FILE__) * "/sim3data.jld") +if !isfile(joinpath(input_dir), "sim3data.jld") centers = randpoints(Ns, dmin, width, height, points) else - import JLD; JLD.@load(dirname(@__FILE__) * "/sim3data.jld", centers) + JLD.@load(joinpath(input_dir, "sim3data.jld"), centers) end N,errN = (934, 9.97040926753751e-7) # @@ -35,26 +34,25 @@ fmm_options = FMMoptions(true, acc = 6, nx = 9, method="pre") divideSpace(centers, fmm_options; drawGroups = false) -draw_fig = false +draw_fig = true # verify and draw begin - assert(verify_min_distance(shapes, centers, ids, points)) + @assert verify_min_distance(shapes, centers, ids, points) if draw_fig figure() #draw shapes and points draw_shapes(shapes, ids, centers, φs) plot(points[:,1], points[:,2], "r*") tight_layout() - ax = gca() - ax[:set_aspect]("equal", adjustable = "box") + axis("equal") end end -plot_border = shapes[1].R*[-1;1;-1;1] + [0.0; width; 0.0; height] +border = shapes[1].R*[-1;1;-1;1] + [0.0; width; 0.0; height] optim_options = Optim.Options(f_tol = 1e-6, - iterations = 100, + iterations = 150, store_trace = true, extended_trace = false, show_trace = true, @@ -62,50 +60,34 @@ optim_options = Optim.Options(f_tol = 1e-6, optim_method = Optim.BFGS(;linesearch = LineSearches.BackTracking()) -tic() -test_max = optimize_φ(φs, points, P, PlaneWave(θ_i), k0, kin, shapes, - centers, ids, fmm_options, optim_options, optim_method; minimize = false) -optim_time = toq() - -# %% +#precompile +test_max = optimize_φ(φs, points, P, ui, k0, kin, shapes, centers, ids, + fmm_options, optim_options, optim_method; minimize = false) +optim_time = @elapsed begin + test_max = optimize_φ(φs, points, P, ui, k0, kin, shapes, centers, ids, + fmm_options, optim_options, optim_method; minimize = false) +end sp_before = ScatteringProblem(shapes, ids, centers, φs) -plot_near_field(k0, kin, P, sp_before, PlaneWave(θ_i), - x_points = 600, y_points = 200, border = plot_border); +Ez1 = plot_near_field(k0, kin, P, sp_before, ui, + x_points = 600, y_points = 200, border = border, method="recurrence") colorbar() clim([0;5]) -plot_near_field_pgf(output_dir * "/opt_phi_before.tex", k0, kin, P, - sp_before, PlaneWave(θ_i); opt = fmm_options, x_points = 150, y_points = 50, - border = plot_border, downsample = 10, include_preamble = true) sp_after = ScatteringProblem(shapes, ids, centers, test_max.minimizer) -plot_near_field(k0, kin, P, sp_after, PlaneWave(θ_i), - x_points = 600, y_points = 200, border = plot_border) +Ez2 = plot_near_field(k0, kin, P, sp_after, ui, + x_points = 600, y_points = 200, border = border, method="recurrence") colorbar() clim([0;5]) -plot_near_field_pgf(output_dir * "/opt_phi_after.tex", k0, kin, P, - sp_after, PlaneWave(θ_i); opt = fmm_options, x_points = 150, y_points = 50, - border = plot_border, downsample = 10, include_preamble = true) inner_iters = length(test_max.trace) fobj = -[test_max.trace[i].value for i=1:inner_iters] gobj = [test_max.trace[i].g_norm for i=1:inner_iters] figure() -plot(0:inner_iters-1, fobj) -plot(0:inner_iters-1, gobj) - -pgf.@pgf begin - fobj_plot = pgf.Plot({blue, thick, no_markers}, - pgf.Coordinates(0:inner_iters-1, fobj)) - gobj_plot = pgf.Plot({red, thick, dashed, no_markers}, - pgf.Coordinates(0:inner_iters-1, gobj)) - ax = pgf.Axis({ width = "6cm", - xlabel = "Iterations", - legend_pos = "north east", - legend_style = "font = \\footnotesize"}, - fobj_plot, gobj_plot) - push!(ax, pgf.Legend(["\$f_{\\mathrm{obj}}\$"; - "\$\\|\\mathbf{g}_{\\mathrm{obj}}\\|\$"])) -end -pgf.save(output_dir * "/opt_phi_conv.tex", ax ,include_preamble = false) +plot(0:inner_iters-1, fobj, label="\$f_{\\mathrm{obj}}\$") +plot(0:inner_iters-1, gobj, label="\$\\Vert\\mathbf{g}_{\\mathrm{obj}}\\Vert\$") +legend(loc="best") +xlabel("Iterations") + +JLD.@save joinpath(output_dir, "angle_optim2.jld") diff --git a/examples/sim3data.jld b/examples/sim3data.jld deleted file mode 100644 index 1ab9dd0..0000000 Binary files a/examples/sim3data.jld and /dev/null differ diff --git a/src/PS_types.jl b/src/PS_types.jl index 7706652..f9b2955 100644 --- a/src/PS_types.jl +++ b/src/PS_types.jl @@ -92,7 +92,7 @@ and the following keyword arguments dictate its behavior: - `dx::Real`: group height/width (alternative division) - `acc::Integer`: accuracy digits for translation truncation, and also for gmres if `tol` is not given - `tol::Real`: gmres tolerance -- `method::String`: method used: for now can be "pre" or "pre2". Mainly used for development. +- `method::String`: method used: for now can be "pre". Mainly used for development. """ mutable struct FMMoptions FMM::Bool #Is FMM used? @@ -100,7 +100,7 @@ mutable struct FMMoptions dx::Real #group height/width (alternative division) acc::Integer #accuracy digits for translation truncation, and also for gmres if tol is not given tol::Real #gmres tolerance - method::String #method used: can be "pre" or "pre2" + method::String #method used: can be "pre" # symmetric::Bool #are agg = disagg points, and thus Disagg[k] = Agg^*[k] ∀k? #empty contructor - for not using FMM @@ -121,8 +121,8 @@ mutable struct FMMoptions error("FMMoptions: accuracy digits must be in [1,16]") tol < 0.0 && error("FMMoptions: gmres tolerance must be greater than 0.0") - !in(method, ("pre","pre2")) && - error("FMMoptions: method must be \"pre\" or \"pre2\"") + !in(method, ["pre"]) && + error("FMMoptions: method must be \"pre\" or ...") tol == 0.0 && (tol = 10^(-Float64(acc))) return new(true, nx, dx, acc, tol, method) end @@ -147,10 +147,10 @@ mutable struct OptimBuffer rhs_grad::Vector{Complex{Float64}} OptimBuffer(Ns::Integer, P::Integer, Npoints::Integer, J = Ns) = - new(Array{Complex{Float64}}(Ns*(2*P+1)), - Array{Complex{Float64}}(Npoints), - Array{Complex{Float64}}(Ns*(2*P+1), J), - Array{Complex{Float64}}(Ns*(2*P+1))) + new(Array{Complex{Float64}}(undef, Ns*(2*P+1)), + Array{Complex{Float64}}(undef, Npoints), + Array{Complex{Float64}}(undef, Ns*(2*P+1), J), + Array{Complex{Float64}}(undef, Ns*(2*P+1))) end """ @@ -214,7 +214,7 @@ start and end points of the source, and `σ` contains the potential density. """ function CurrentSource(x1, y1, x2, y2, σ) len = sqrt((x2 - x1)^2 + (y2 - y1)^2) - t = linspace(0, 1, length(σ)) + t = range(0, stop=1, length=length(σ)) p = [(x2 - x1)*t + x1 (y2 - y1)*t + y1] CurrentSource(p, σ, len) end diff --git a/src/ParticleScattering.jl b/src/ParticleScattering.jl index 5b271e3..5815746 100644 --- a/src/ParticleScattering.jl +++ b/src/ParticleScattering.jl @@ -7,11 +7,15 @@ for various design problems. module ParticleScattering #Core functionality -using SpecialFunctions, IterativeSolvers, LinearMaps, Optim +using SpecialFunctions, IterativeSolvers, LinearMaps, Optim, SparseArrays, LinearAlgebra +import Statistics: mean import LineSearches -#For plotting with PyPlot +#For plotting with PyPlot - actual import is done at runtime using PyPlot, PyCall -@pyimport matplotlib.patches as patch #circles, polygons +const patch = PyNULL() +function __init__() + copy!(patch, pyimport("matplotlib.patches")) #circles, polygons +end include("PS_types.jl") include("shapes.jl") diff --git a/src/fmm_main.jl b/src/fmm_main.jl index 325733a..6ba804b 100644 --- a/src/fmm_main.jl +++ b/src/fmm_main.jl @@ -11,88 +11,78 @@ or inner cylindrical coefficients (in case of circular). `plot_res` controls plotting of the residual. Inner coefficients are calculated only if `get_inner` is true, and timing is printed if `verbose` is true. """ -#TODO: return beta,inner,history function solve_particle_scattering_FMM(k0, kin, P, sp::ScatteringProblem, u::Einc, opt::FMMoptions; plot_res = false, get_inner = true, verbose = true) - assert(opt.FMM) + @assert opt.FMM "opt.FMM is disabled, set to true" shapes = sp.shapes; ids = sp.ids; centers = sp.centers; φs = sp.φs Ns = size(sp) groups, boxSize = divideSpace(centers, opt) P2, Q = FMMtruncation(opt.acc, boxSize, k0) verbose && println("FMM solution timing:") - tic() - mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri=true) - dt0 = toq() - - tic() - scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids) - dt1 = toq() + dt0 = @elapsed begin + mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri=true) + end - tic() - #construct rhs - rhs = u2α(k0, u, centers, P) - for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) - #see if there is a faster alternative - if φs[ic] == 0.0 - rhs[rng] = scatteringMatrices[ids[ic]]*rhs[rng] - else - #rotate without matrix - rotateMultipole!(view(rhs,rng),-φs[ic],P) - rhs[rng] = scatteringMatrices[ids[ic]]*rhs[rng] - rotateMultipole!(view(rhs,rng),φs[ic],P) - end + dt1 = @elapsed begin + scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids) end - dt2 = toq() + + dt2 = @elapsed begin + #construct rhs + rhs = u2α(k0, u, centers, P) + for ic = 1:Ns + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) + #see if there is a faster alternative + if φs[ic] == 0.0 + rhs[rng] = scatteringMatrices[ids[ic]]*rhs[rng] + else + #rotate without matrix + rotateMultipole!(view(rhs,rng),-φs[ic],P) + rhs[rng] = scatteringMatrices[ids[ic]]*rhs[rng] + rotateMultipole!(view(rhs,rng),φs[ic],P) + end + end + end pre_agg_buffer = zeros(Complex{Float64},Q,length(groups)) - trans_buffer = Array{Complex{Float64}}(Q) - tic() - if opt.method == "pre" - MVP = LinearMap{eltype(rhs)}((output_, x_) -> FMM_mainMVP_pre!(output_, + trans_buffer = Array{Complex{Float64}}(undef, Q) + dt3 = @elapsed begin + MVP = LinearMap{eltype(rhs)}((output_, x_) -> FMM_mainMVP_pre!(output_, x_, scatteringMatrices, φs, ids, P, mFMM, pre_agg_buffer, trans_buffer), Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - elseif opt.method == "pre2" - MVP = LinearMap{eltype(rhs)}((output_, x_) -> FMM_mainMVP_pre2!(output_, - x_, scatteringMatrices, φs, ids, P, mFMM, - pre_agg_buffer, trans_buffer), Ns*(2*P+1), - Ns*(2*P+1), ismutating = true) + result = gmres(MVP, rhs, restart = Ns*(2*P+1), tol = opt.tol, log = true) #no restart, preconditioning end - result = gmres(MVP, rhs, restart = Ns*(2*P+1), tol = opt.tol, log = true) #no restart, preconditioning - dt3 = toq() if get_inner #recover full incoming expansion - in sigma_mu terms for parametrized shape, #in multipole expansion for circle - tic() - - #find LU factorization once for each shape - scatteringLU = [lufact(scatteringMatrices[i]) for i = 1:length(shapes)] - - inner = Array{Vector{Complex{Float64}}}(Ns) - α_c = Array{Complex{Float64}}(2*P+1) - for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) - if typeof(shapes[ids[ic]]) == ShapeParams - if φs[ic] == 0.0 - α_c[:] = scatteringLU[ids[ic]]\result[1][rng] - else - rotateMultipole!(α_c,view(result[1],rng),-φs[ic],P) - A_ldiv_B!(scatteringLU[ids[ic]],α_c) - end - inner[ic] = innerExpansions[ids[ic]]*α_c - else - inner[ic] = innerExpansions[ids[ic]]*result[1][rng] - end - end - dt4 = toq() + dt4 = @elapsed begin + #find LU factorization once for each shape + scatteringLU = [lu(scatteringMatrices[i]) for i = 1:length(shapes)] + inner = Array{Vector{Complex{Float64}}}(undef, Ns) + α_c = Array{Complex{Float64}}(undef, 2*P+1) + for ic = 1:Ns + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) + if typeof(shapes[ids[ic]]) == ShapeParams + if φs[ic] == 0.0 + α_c[:] = scatteringLU[ids[ic]]\result[1][rng] + else + rotateMultipole!(α_c, view(result[1],rng), -φs[ic], P) + ldiv!(scatteringLU[ids[ic]], α_c) + end + inner[ic] = innerExpansions[ids[ic]]*α_c + else + inner[ic] = innerExpansions[ids[ic]]*result[1][rng] + end + end + end end if plot_res residual_ = MVP*result[1] - rhs println("last residual is empirically: $(norm(residual_)), but result[2].residuals[end] = $(result[2].residuals[end])") figure() - semilogy(result[2].residuals.') + semilogy(transpose(result[2].residuals)) end if verbose println("FMM matrix construction: $dt0 s") diff --git a/src/fmm_matrices.jl b/src/fmm_matrices.jl index 2a66f4d..666ef95 100644 --- a/src/fmm_matrices.jl +++ b/src/fmm_matrices.jl @@ -9,8 +9,8 @@ end function divideSpace(centers::Array{Float64,2}, options; drawGroups = false) size(centers,1) < 2 && error("divideSpace: need at least 2 points") # using findfirst everywhere is slower but more robust to floating errors - (x_min,y_min) = minimum(centers,1) - (x_max,y_max) = maximum(centers,1) + (x_min,y_min) = minimum(centers, dims=1) + (x_max,y_max) = maximum(centers, dims=1) epss = eps(maximum(abs.([x_min;x_max;y_min;y_max]))) x_max += epss y_max += epss @@ -37,7 +37,7 @@ function divideSpace(centers::Array{Float64,2}, options; drawGroups = false) leftover_y = max(a*ny - (y_max-y_min),0.0) y_max += 0.5*leftover_y y_min -= 0.5*leftover_y - FMMgroups = Array{FMMgroup}(nx,ny) + FMMgroups = Array{FMMgroup}(undef,nx,ny) for ix = 1:nx, iy = 1:ny center_x = x_min + (ix-0.5)*a center_y = y_min + (iy-0.5)*a @@ -70,7 +70,7 @@ function divideSpace(centers::Array{Float64,2}, options; drawGroups = false) ylim([y_min - a;y_max + a]) tight_layout() ax = gca() - ax[:set_aspect]("equal", adjustable = "box") + ax.set_aspect("equal", adjustable = "box") end return FMMgroups,a end @@ -82,7 +82,7 @@ end # for i = 1:size(centers,1) # d[1] = centers[i,1] - box_center[1] # d[2] = centers[i,2] - box_center[2] -# A[:,(i-1)*(2*P+1) + (1:2*P+1)] = [exp(-1.0im*(k*(cos(t[q])*d[1] + sin(t[q])*d[2]) +# A[:,(i-1)*(2*P+1) .+ (1:2*P+1)] = [exp(-1.0im*(k*(cos(t[q])*d[1] + sin(t[q])*d[2]) # + n*(pi/2-t[q]))) for q=1:Q, n=-P:P] # end # end @@ -91,12 +91,12 @@ function FMMaggregationMatrix(k, centers, box_center, t, P) #P is length of first multipole expansion, *not* FMM Q = length(t) Ns = size(centers,1) - A = Array{Complex{Float64}}(Q, Ns*(2*P+1)) - d = Array{Float64}(2) #silly but faster + A = Array{Complex{Float64}}(undef, Q, Ns*(2*P+1)) + d = Array{Float64}(undef, 2) #silly but faster for i = 1:Ns d[1] = centers[i,1] - box_center[1] d[2] = centers[i,2] - box_center[2] - A[:,(i-1)*(2*P+1) + (1:2*P+1)] = [exp(-1.0im*(k*(cos(t[q])*d[1] + sin(t[q])*d[2]) + A[:,(i-1)*(2*P+1) .+ (1:2*P+1)] = [exp(-1.0im*(k*(cos(t[q])*d[1] + sin(t[q])*d[2]) + n*(pi/2-t[q]))) for q=1:Q, n=-P:P] end return A @@ -119,9 +119,9 @@ function FMMtranslationMatrix(k, x, t, P2) Q = length(t) T = zeros(Complex{Float64},Q) nx = sqrt(sum(abs2,x)) - tx = Float64[-atan2(-x[2],-x[1]) - pi/2 + t[q] for q=1:Q] + tx = Float64[-atan(-x[2],-x[1]) - pi/2 + t[q] for q=1:Q] bess = besselh.(-P2:P2, 1, k*nx) - T[:] = 0 + T[:] .= 0 for i = -P2:P2, j = 1:Q T[j] += bess[i+P2+1]*exp(1.0im*i*tx[j])/Q end @@ -133,35 +133,37 @@ function FMMnearMatrix(k, P, groups, centers, boxSize, num) Ns = size(centers,1) G = length(groups) W = 2*P+1 - Is = Array{Int}(num*W^2) - Js = Array{Int}(num*W^2) - Zs = Array{Complex{Float64}}(num*W^2) - bess = Array{Complex{Float64}}(2*P+1) + Is = Array{Int}(undef, num*W^2) + Js = Array{Int}(undef, num*W^2) + Zs = Array{Complex{Float64}}(undef, num*W^2) + bess = Array{Complex{Float64}}(undef, 2*P+1) mindist2 = 3*boxSize^2 #anywhere between 2 and 4 ind = 0 - for ig1 = 1:G, ig2 = 1:G - x = groups[ig1].center - groups[ig2].center - sum(abs2,x) > mindist2 && continue - #for each enclosed scatterer, build translation matrix - for ic1 = 1:length(groups[ig1].point_ids) - for ic2 = 1:length(groups[ig2].point_ids) - ig1 == ig2 && (ic1 == ic2 && continue) #no self-interactions - d = centers[groups[ig1].point_ids[ic1],:] - centers[groups[ig2].point_ids[ic2],:] - nd = sqrt(sum(abs2,d)) - td = atan2(d[2],d[1]) - bess[:] = besselh.(0:2*P,1,k*nd) - for ix = 1:2*P #lower diagonals - rng = ix+1:1+W:W^2-W*ix - Zs[ind + rng] = -exp(-1im*td*ix)*(-1)^(ix)*bess[ix+1] - end - for ix = 0:2*P #central and upper diagonals - rng = ix*W+1:1+W:W^2-ix - Zs[ind + rng] = -exp(1im*td*ix)*bess[ix+1] - end - for ij = 1:W - Is[ind + (1:W)] = collect((groups[ig1].point_ids[ic1]-1)*W + (1:W)) - Js[ind + (1:W)] = (groups[ig2].point_ids[ic2]-1)*W + ij - ind += W + for ig1 = 1:G + for ig2 = 1:G + x = groups[ig1].center - groups[ig2].center + sum(abs2,x) > mindist2 && continue + #for each enclosed scatterer, build translation matrix + for ic1 = 1:length(groups[ig1].point_ids) + for ic2 = 1:length(groups[ig2].point_ids) + ig1 == ig2 && (ic1 == ic2 && continue) #no self-interactions + d = centers[groups[ig1].point_ids[ic1],:] - centers[groups[ig2].point_ids[ic2],:] + nd = sqrt(sum(abs2,d)) + td = atan(d[2],d[1]) + bess[:] = besselh.(0:2*P,1,k*nd) + for ix = 1:2*P #lower diagonals + rng = ix+1:1+W:W^2-W*ix + Zs[ind .+ rng] .= -exp(-1im*td*ix)*(-1)^(ix)*bess[ix+1] + end + for ix = 0:2*P #central and upper diagonals + rng = ix*W+1:1+W:W^2-ix + Zs[ind .+ rng] .= -exp(1im*td*ix)*bess[ix+1] + end + for ij = 1:W + Is[ind .+ (1:W)] = collect((groups[ig1].point_ids[ic1]-1)*W .+ (1:W)) + Js[ind .+ (1:W)] .= (groups[ig2].point_ids[ic2]-1)*W + ij + ind += W + end end end end @@ -175,13 +177,13 @@ function FMMnearMatrix_upperTri(k, P, groups, centers, boxSize, num) Ns = size(centers,1) G = length(groups) W = 2*P+1 - Is = Array{Int}(num*W^2) - Js = Array{Int}(num*W^2) - Zs = Array{Complex{Float64}}(num*W^2) + Is = Array{Int}(undef, num*W^2) + Js = Array{Int}(undef, num*W^2) + Zs = Array{Complex{Float64}}(undef, num*W^2) mindist2 = 3*boxSize^2 #anywhere between 2 and 4 ind = 0 - d = Array{Float64}(2) - bess = Array{Complex{Float64}}(2*P+1) + d = Array{Float64}(undef, 2) + bess = Array{Complex{Float64}}(undef, 2*P+1) #first, between group and itself for ig1 = 1:G for ic1 = 1:length(groups[ig1].point_ids) @@ -190,74 +192,76 @@ function FMMnearMatrix_upperTri(k, P, groups, centers, boxSize, num) d[2] = centers[groups[ig1].point_ids[ic1],2] - centers[groups[ig1].point_ids[ic2],2] nd = sqrt(sum(abs2,d)) #first for ic2->ic1, then ic1->ic2 - td1 = atan2(d[2],d[1]) - td2 = atan2(-d[2],-d[1]) #adding pi instead leads to error + td1 = atan(d[2],d[1]) + td2 = atan(-d[2],-d[1]) #adding pi instead leads to error bess[:] = besselh.(0:2*P,1,k*nd) ind2 = ind + W^2 for ix = 1:2*P #lower diagonals rng = ix+1:1+W:W^2-W*ix - Zs[ind + rng] = -exp(-1.0im*td1*ix)*(-1)^(ix)*bess[ix+1] - Zs[ind2 + rng] = -exp(-1.0im*td2*ix)*(-1)^(ix)*bess[ix+1] + Zs[ind .+ rng] .= -exp(-1.0im*td1*ix)*(-1)^(ix)*bess[ix+1] + Zs[ind2 .+ rng] .= -exp(-1.0im*td2*ix)*(-1)^(ix)*bess[ix+1] end for ix = 0:2*P #central and upper diagonals rng = ix*W+1:1+W:W^2-ix - Zs[ind + rng] = -exp(1.0im*td1*ix)*bess[ix+1] - Zs[ind2 + rng] = -exp(1.0im*td2*ix)*bess[ix+1] + Zs[ind .+ rng] .= -exp(1.0im*td1*ix)*bess[ix+1] + Zs[ind2 .+ rng] .= -exp(1.0im*td2*ix)*bess[ix+1] end #this way we avoid allocating vectors for ij = 1:W - rng = (1:W) + (ij-1)*W - Js[ind + rng] = (groups[ig1].point_ids[ic2]-1)*W + ij - Js[ind2 + rng] = (groups[ig1].point_ids[ic1]-1)*W + ij + rng = (1:W) .+ (ij-1)*W + Js[ind .+ rng] .= (groups[ig1].point_ids[ic2]-1)*W + ij + Js[ind2 .+ rng] .= (groups[ig1].point_ids[ic1]-1)*W + ij end for ij = 1:W - rng = (1:W:W^2-W+1) + (ij-1) - Is[ind + rng] = (groups[ig1].point_ids[ic1]-1)*W + ij - Is[ind2 + rng] = (groups[ig1].point_ids[ic2]-1)*W + ij + rng = (1:W:W^2-W+1) .+ (ij-1) + Is[ind .+ rng] .= (groups[ig1].point_ids[ic1]-1)*W + ij + Is[ind2 .+ rng] .= (groups[ig1].point_ids[ic2]-1)*W + ij end ind += 2*W^2 #skip over ind2 data end end end #now, between groups - for ig1 = 1:G, ig2 = ig1+1:G - x = groups[ig1].center - groups[ig2].center - sum(abs2,x) > mindist2 && continue - #for each enclosed scatterer, build translation matrix - for ic1 = 1:length(groups[ig1].point_ids) - for ic2 = 1:length(groups[ig2].point_ids) - d[1] = centers[groups[ig1].point_ids[ic1],1] - centers[groups[ig2].point_ids[ic2],1] - d[2] = centers[groups[ig1].point_ids[ic1],2] - centers[groups[ig2].point_ids[ic2],2] - nd = sqrt(sum(abs2,d)) - #first for ig2->ig1, then ig1->ig2 - td1 = atan2(d[2],d[1]) - td2 = atan2(-d[2],-d[1]) #adding pi instead leads to error - bess[:] = besselh.(0:2*P,1,k*nd) - ind2 = ind + W^2 + for ig1 = 1:G + for ig2 = ig1+1:G + x = groups[ig1].center - groups[ig2].center + sum(abs2,x) > mindist2 && continue + #for each enclosed scatterer, build translation matrix + for ic1 = 1:length(groups[ig1].point_ids) + for ic2 = 1:length(groups[ig2].point_ids) + d[1] = centers[groups[ig1].point_ids[ic1],1] - centers[groups[ig2].point_ids[ic2],1] + d[2] = centers[groups[ig1].point_ids[ic1],2] - centers[groups[ig2].point_ids[ic2],2] + nd = sqrt(sum(abs2,d)) + #first for ig2->ig1, then ig1->ig2 + td1 = atan(d[2],d[1]) + td2 = atan(-d[2],-d[1]) #adding pi instead leads to error + bess[:] = besselh.(0:2*P,1,k*nd) + ind2 = ind + W^2 - for ix = 1:2*P #lower diagonals - rng = ix+1:1+W:W^2-W*ix - Zs[ind + rng] = -exp(-1.0im*td1*ix)*(-1)^(ix)*bess[ix+1] - Zs[ind2 + rng] = -exp(-1.0im*td2*ix)*(-1)^(ix)*bess[ix+1] - end - for ix = 0:2*P #central and upper diagonals - rng = ix*W+1:1+W:W^2-ix - Zs[ind + rng] = -exp(1.0im*td1*ix)*bess[ix+1] - Zs[ind2 + rng] = -exp(1.0im*td2*ix)*bess[ix+1] - end - #this way we avoid allocating vectors - for ij = 1:W - rng = (1:W) + (ij-1)*W - Js[ind + rng] = (groups[ig2].point_ids[ic2]-1)*W + ij - Js[ind2 + rng] = (groups[ig1].point_ids[ic1]-1)*W + ij - end - for ij = 1:W - rng = (1:W:W^2-W+1) + (ij-1) - Is[ind + rng] = (groups[ig1].point_ids[ic1]-1)*W + ij - Is[ind2 + rng] = (groups[ig2].point_ids[ic2]-1)*W + ij + for ix = 1:2*P #lower diagonals + rng = ix+1:1+W:W^2-W*ix + Zs[ind .+ rng] .= -exp(-1.0im*td1*ix)*(-1)^(ix)*bess[ix+1] + Zs[ind2 .+ rng] .= -exp(-1.0im*td2*ix)*(-1)^(ix)*bess[ix+1] + end + for ix = 0:2*P #central and upper diagonals + rng = ix*W+1:1+W:W^2-ix + Zs[ind .+ rng] .= -exp(1.0im*td1*ix)*bess[ix+1] + Zs[ind2 .+ rng] .= -exp(1.0im*td2*ix)*bess[ix+1] + end + #this way we avoid allocating vectors + for ij = 1:W + rng = (1:W) .+ (ij-1)*W + Js[ind .+ rng] .= (groups[ig2].point_ids[ic2]-1)*W + ij + Js[ind2 .+ rng] .= (groups[ig1].point_ids[ic1]-1)*W + ij + end + for ij = 1:W + rng = (1:W:W^2-W+1) .+ (ij-1) + Is[ind .+ rng] .= (groups[ig1].point_ids[ic1]-1)*W + ij + Is[ind2 .+ rng] .= (groups[ig2].point_ids[ic2]-1)*W + ij + end + ind += 2*W^2 #skip over ind2 data end - ind += 2*W^2 #skip over ind2 data end end end @@ -278,7 +282,7 @@ function FMMbuildMatrices(k, P, P2, Q, groups, centers, boxSize; tri = true) Trans = Vector{Complex{Float64}}[] #For now, the translation from j to i is at Trans[(i-1)*G + j] mindist2 = 3*boxSize^2 #anywhere between 2 and 4 - d = Array{Float64}(2) + d = Array{Float64}(undef, 2) num = 0 for ig1 = 1:G for ig2 = 1:G diff --git a/src/fmm_mvp.jl b/src/fmm_mvp.jl index 8e74529..59d44b1 100644 --- a/src/fmm_mvp.jl +++ b/src/fmm_mvp.jl @@ -3,7 +3,7 @@ function FMM_mainMVP_pre!(output, beta, scatteringMatrices, φs::Vector{Float64} #calculate matrix-vector product - devectorized with pre-preconditioning @inbounds begin - A_mul_B!(output, mFMM.Znear, beta) + mul!(output, mFMM.Znear, beta) G = length(mFMM.groups) fill!(pre_agg,0.0) #first preagg all @@ -45,90 +45,20 @@ function FMM_mainMVP_pre!(output, beta, scatteringMatrices, φs::Vector{Float64} #multiply by S to produce -S #temp can be moved outward, but for now this preallocation prevents most #dynamic mem alloc by this MVP - temp = Array{Complex{Float64}}(2*P+1) + temp = Array{Complex{Float64}}(undef, 2*P+1) for ic = 1:length(ids) - rng = (ic-1)*(2*P+1) + (1:2*P+1) - #copy values to avoid dynamic allocation - # copy!(temp,1:2*P+1,output,rng) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) + #copy values to avoid dynamic allocation - copyto!(temp,1:2*P+1,output,rng) doesn't work because output is a subarray for ip = 1:2*P+1 temp[ip] = output[(ic-1)*(2*P+1) + ip] end v = view(output,rng) if φs[ic] == 0.0 - A_mul_B!(v, scatteringMatrices[ids[ic]], temp) + mul!(v, scatteringMatrices[ids[ic]], temp) else #rotate without matrix rotateMultipole!(temp,-φs[ic],P) - A_mul_B!(v, scatteringMatrices[ids[ic]], temp) - rotateMultipole!(v,φs[ic],P) - end - end - #add identity matrix (Ix) - output .+= beta - end #inbounds - return output -end - -function FMM_mainMVP_pre2!(output, beta, scatteringMatrices, φs::Vector{Float64}, ids::Vector{Int}, P, mFMM, pre_agg, translated_sum) - #calculate matrix-vector product - partially vectorized with pre-preconditioning - - @inbounds begin - A_mul_B!(output, mFMM.Znear, beta) - G = length(mFMM.groups) - fill!(pre_agg,0.0) - #first preagg all - for ig2 = 1:G - for is = 1:mFMM.groups[ig2].size - indices = (mFMM.groups[ig2].point_ids[is]-1)*(2*P+1) - for ii = 1:2*P+1 - for iq = 1:mFMM.Q - pre_agg[iq,ig2] += mFMM.Agg[ig2][iq,(is-1)*(2*P+1) + ii]*beta[indices + ii] - end - end - end - end - - for ig1 = 1:G - #translate plane waves from ig2 to ig1 - fill!(translated_sum,0.0) - for ig2 = 1:G - if isempty(mFMM.Trans[(ig1-1)*G+ig2])#FMMnear(ig1,ig2) - either compute distance or check if Trans is empty - continue - else - for iQ = 1:mFMM.Q - translated_sum[iQ] -= mFMM.Trans[(ig1-1)*G + ig2][iQ]*pre_agg[iQ,ig2] - end - #minus sign because the real equation is (I-ST)x=b - end - end - #disaggregate from ig1 center to ig1's scatterers, producing -Tx - for is = 1:mFMM.groups[ig1].size - for ip = 1:2*P+1 - disagged = 0.0im - for iq = 1:mFMM.Q - disagged += conj(mFMM.Agg[ig1][iq,(is-1)*(2*P+1) + ip])*translated_sum[iq] - end - output[(mFMM.groups[ig1].point_ids[is]-1)*(2*P+1) + ip] += disagged - end - end - end - #multiply by S to produce -ST - #temp can be moved outward, but for now this preallocation prevents most - #dynamic mem alloc by this MVP - temp = Array{Complex{Float64}}(2*P+1) - for ic = 1:length(ids) - rng = (ic-1)*(2*P+1) + (1:2*P+1) - #copy values to avoid dynamic allocation - for ip = 1:2*P+1 - temp[ip] = output[(ic-1)*(2*P+1) + ip] - end - v = view(output,rng) - if φs[ic] == 0.0 - A_mul_B!(v, scatteringMatrices[ids[ic]], temp) - else - #rotate without matrix - rotateMultipole!(temp,-φs[ic],P) - A_mul_B!(v, scatteringMatrices[ids[ic]], temp) + mul!(v, scatteringMatrices[ids[ic]], temp) rotateMultipole!(v,φs[ic],P) end end @@ -146,22 +76,22 @@ function FMM_mainMVP_transpose!(output, beta, scatteringMatrices, φs::Vector{Fl #Xβ storage can be moved outward, but for now this preallocation prevents most #dynamic mem alloc by this MVP - Xβ = Array{Complex{Float64}}(length(beta)) - temp = Array{Complex{Float64}}(2*P+1) + Xβ = Array{Complex{Float64}}(undef, length(beta)) + temp = Array{Complex{Float64}}(undef, 2*P+1) for ic = 1:length(ids) - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) v = view(Xβ, rng) if φs[ic] == 0.0 - At_mul_B!(v, scatteringMatrices[ids[ic]], view(beta, rng)) + mul!(v, transpose(scatteringMatrices[ids[ic]]), view(beta, rng)) else #rotate without matrix - transposed rotateMultipole!(temp, view(beta, rng), φs[ic], P) - At_mul_B!(v, scatteringMatrices[ids[ic]], temp) + mul!(v, transpose(scatteringMatrices[ids[ic]]), temp) rotateMultipole!(v, -φs[ic], P) end end - At_mul_B!(output, mFMM.Znear, Xβ) + mul!(output, transpose(mFMM.Znear), Xβ) G = length(mFMM.groups) fill!(pre_agg,0.0) #first preagg all diff --git a/src/incident.jl b/src/incident.jl index a415973..efc9ea1 100644 --- a/src/incident.jl +++ b/src/incident.jl @@ -3,7 +3,7 @@ const eta0 = 4π*299792458e-7 #plane wave function u2α(k, u::PlaneWave, centers::Array{T,2}, P) where T <: Real #build incoming coefficients for plane wave incident field - α = Array{Complex{Float64}}(size(centers,1)*(2P + 1)) + α = Array{Complex{Float64}}(undef, size(centers,1)*(2P + 1)) for ic = 1:size(centers,1) phase = exp(1.0im*k*(cos(u.θi)*centers[ic,1] + sin(u.θi)*centers[ic,2])) #phase shift added to move cylinder coords for p = -P:P @@ -28,10 +28,10 @@ end #line (filament) source function u2α(k, u::LineSource, centers::Array{T,2}, P) where T <: Real - α = Array{Complex{Float64}}(size(centers,1)*(2P + 1)) + α = Array{Complex{Float64}}(undef, size(centers,1)*(2P + 1)) for ic = 1:size(centers,1) R = hypot(centers[ic,1] - u.x, centers[ic,2] - u.y) - θ = atan2(centers[ic,2] - u.y, centers[ic,1] - u.x) + θ = atan(centers[ic,2] - u.y, centers[ic,1] - u.x) α[(ic-1)*(2P+1) + 0 + P + 1] = 0.25im*besselh(0, 1, k*R) for p = 1:P α[(ic-1)*(2P+1) - p + P + 1] = 0.25im*exp(1im*p*θ)*besselh(p, 1, k*R) @@ -42,9 +42,9 @@ function u2α(k, u::LineSource, centers::Array{T,2}, P) where T <: Real end function uinc(k, points::Array{T,2}, u::LineSource) where T <: Real - r = hypot.(u.x - points[:,1], u.y - points[:,2]) + r = hypot.(u.x .- points[:,1], u.y .- points[:,2]) if any(r .== 0) - warn("uinc: encountered singularity in incident field, returned NaN") + @warn("uinc: encountered singularity in incident field, returned NaN") r[r.==0] = NaN end 0.25im*besselh.(0, 1, k*r) @@ -53,7 +53,7 @@ end function uinc(k, points::Array{T,1}, u::LineSource) where T <: Real r = hypot(u.x - points[1], u.y - points[2]) if r == 0 - warn("uinc: encountered singularity in incident field, returned NaN") + @warn("uinc: encountered singularity in incident field, returned NaN") return NaN end 0.25im*besselh(0, 1, k*r) @@ -66,7 +66,7 @@ function u2α(k, u::CurrentSource, centers::Array{T,2}, P) where T <: Real for ic = 1:size(centers,1) for is = 1:length(u.σ) R = hypot(centers[ic,1] - u.p[is,1], centers[ic,2] - u.p[is,2]) - θ = atan2(centers[ic,2] - u.p[is,2], centers[ic,1] - u.p[is,1]) + θ = atan(centers[ic,2] - u.p[is,2], centers[ic,1] - u.p[is,1]) α[(ic-1)*(2P+1) + P + 1] += c*u.σ[is]*besselh(0, 1, k*R) for p = 1:P @@ -80,12 +80,12 @@ function u2α(k, u::CurrentSource, centers::Array{T,2}, P) where T <: Real end function uinc(k, points::Array{T,2}, u::CurrentSource) where T <: Real - res = Array{Complex{Float64}}(size(points, 1)) - r = Array{Float64}(size(u.p, 1)) + res = Array{Complex{Float64}}(undef, size(points, 1)) + r = Array{Float64}(undef, size(u.p, 1)) for i = 1:size(points, 1) r .= hypot.(u.p[:,1] - points[i,1], u.p[:,2] - points[i,2]) if any(r .== 0) #point is on source - warn("uinc: encountered singularity in incident field, returned NaN") + @warn("uinc: encountered singularity in incident field, returned NaN") res[i] = NaN end res[i] = (0.25im*u.len/length(u.σ))*sum(besselh.(0, 1, k*r).*u.σ) @@ -96,7 +96,7 @@ end function uinc(k, points::Array{T,1}, u::CurrentSource) where T <: Real r = hypot.(u.p[:,1] - points[1], u.p[:,2] - points[2]) if any(r .== 0) #point is on source - warn("uinc: encountered singularity in incident field, returned NaN") + @warn("uinc: encountered singularity in incident field, returned NaN") return NaN end (0.25im*u.len/length(u.σ))*sum(besselh.(0, 1, k*r).*u.σ) @@ -119,7 +119,7 @@ function err_α(k, ui::Einc, P, shape, center) end Nt = 10_000 - θ = linspace(0, 2π, Nt)[1:end-1] + θ = range(0, stop=2π, length=Nt)[1:end-1] α = u2α(k, ui, center, P) u_exact = uinc(k, center .+ shape.R*[cos.(θ) sin.(θ)], ui) u_α = sum(α[p + P + 1]*besselj(p, k*shape.R)*exp.(1im*p*θ) for p = -P:P) @@ -136,33 +136,33 @@ function hyinc(k, p, u::PlaneWave) end function hxinc(k, p, u::LineSource) - R = sqrt.((view(p,:,1) - u.x).^2 + (view(p,:,2) - u.y).^2) + R = sqrt.((view(p,:,1) .- u.x).^2 + (view(p,:,2) .- u.y).^2) if any(R .== 0) - warn("Hinc: encountered singularity in incident field, returned NaN") + @warn("Hinc: encountered singularity in incident field, returned NaN") R[R.==0] = NaN end - h = (0.25/eta0./R).*besselh.(1,k*R).*(u.y - view(p,:,2)) + h = (0.25/eta0./R).*besselh.(1,k*R).*(u.y .- view(p,:,2)) end function hyinc(k, p, u::LineSource) - R = sqrt.((view(p,:,1) - u.x).^2 + (view(p,:,2) - u.y).^2) + R = sqrt.((view(p,:,1) .- u.x).^2 + (view(p,:,2) .- u.y).^2) if any(R .== 0) - warn("Hinc: encountered singularity in incident field, returned NaN") + @warn("Hinc: encountered singularity in incident field, returned NaN") R[R.==0] = NaN end - h = (0.25/eta0./R).*besselh.(1,k*R).*(view(p,:,1) - u.x) + h = (0.25/eta0./R).*besselh.(1,k*R).*(view(p,:,1) .- u.x) end #these untested function hxinc(k, points, u::CurrentSource) - res = Array{Complex{Float64}}(size(points, 1)) - r = Array{Float64}(size(u.p, 1)) - y = Array{Float64}(size(u.p, 1)) + res = Array{Complex{Float64}}(undef, size(points, 1)) + r = Array{Float64}(undef, size(u.p, 1)) + y = Array{Float64}(undef, size(u.p, 1)) for i = 1:size(points, 1) - y .= points[i,2] - u.p[:,2] - r .= hypot.(points[i,1] - u.p[:,1], y) + y[:] = points[i,2] - view(u.p,:,2) + r[:] = hypot.(points[i,1] .- view(u.p,:,1), y) if any(r .== 0) #point is on source - warn("Hinc: encountered singularity in incident field, returned NaN") + @warn("Hinc: encountered singularity in incident field, returned NaN") res[i] = NaN end res[i] = (-0.25*u.len/length(u.σ)/eta0)*sum(besselh.(1, 1, k*r).*u.σ.*y./r) @@ -170,14 +170,14 @@ function hxinc(k, points, u::CurrentSource) res end function hyinc(k, points, u::CurrentSource) - res = Array{Complex{Float64}}(size(points, 1)) - r = Array{Float64}(size(u.p, 1)) - x = Array{Float64}(size(u.p, 1)) + res = Array{Complex{Float64}}(undef, size(points, 1)) + r = Array{Float64}(undef, size(u.p, 1)) + x = Array{Float64}(undef, size(u.p, 1)) for i = 1:size(points, 1) - x .= points[i,1] - u.p[:,1] - r .= hypot.(x, points[i,2] - u.p[:,2]) + x[:] = points[i,1] .- view(u.p,:,1) + r[:] = hypot.(x, points[i,2] .- view(u.p,:,2)) if any(r .== 0) #point is on source - warn("Hinc: encountered singularity in incident field, returned NaN") + @warn("Hinc: encountered singularity in incident field, returned NaN") res[i] = NaN end res[i] = (0.25*u.len/length(u.σ)/eta0)*sum(besselh.(1, 1, k*r).*u.σ.*x./r) diff --git a/src/minimum_N_P.jl b/src/minimum_N_P.jl index 5197001..7b5a904 100644 --- a/src/minimum_N_P.jl +++ b/src/minimum_N_P.jl @@ -22,7 +22,7 @@ function minimumN(kout, kin, shape_function; tol = 1e-9, N_points = 10_000, N_st s = shape_function(N_max) #just for radius err_points = [s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)] E_ana = (0.25im*besselh(0, 1, kout*s.R))*ones(Complex{Float64}, N_points) - E_comp = Array{Complex{Float64}}(N_points) + E_comp = Array{Complex{Float64}}(undef, N_points) N = N_start err_start = minimumN_helper(N, kout, kin, shape_function, err_points, E_comp, E_ana) @@ -36,18 +36,18 @@ function minimumN(kout, kin, shape_function; tol = 1e-9, N_points = 10_000, N_st if err > tol if N > N_max - warn("minimumN: Cannot enter binary search mode.") + @warn("minimumN: Cannot enter binary search mode.") else - warn("minimumN: Entering binary search mode in [$N,$N_max].") + @warn("minimumN: Entering binary search mode in [$N,$N_max].") N, err = binarySearch(N_ -> -minimumN_helper(N_, kout, kin, shape_function, err_points, E_comp, E_ana), -tol, N, N_max) end - err > tol && warn("minimumN: Failed to find err(N) < tol.") + err > tol && @warn("minimumN: Failed to find err(N) < tol.") elseif err < tol if N == N_min return N,err end - warn("minimumN: Entering binary search in [$N_min,$N].") + @warn("minimumN: Entering binary search in [$N_min,$N].") N, err = binarySearch(N_ -> -minimumN_helper(N_, kout, kin, shape_function, err_points, E_comp, E_ana), -tol, N_min, N) end @@ -79,12 +79,12 @@ Uses binary search between `P_min` and `P_max`. """ function minimumP(k0, kin, s::ShapeParams; tol = 1e-9, N_points = 10_000, P_min = 1, P_max = 60, dist = 2.0) err_points = [dist*s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)] - warn("""minimumP: Calculating at $dist*R, thus implicitly assuming all scatterers have same radius and we do not care about closer implications.""") + @warn("""minimumP: Calculating at $dist*R, thus implicitly assuming all scatterers have same radius and we do not care about closer implications.""") #compute direct solution for comparison inner = get_potentialPW(k0, kin, s, 0.0) E_quadrature = scatteredfield(inner, k0, s, err_points) - E_multipole = Array{Complex{Float64}}(N_points) + E_multipole = Array{Complex{Float64}}(undef, N_points) P, errP = binarySearch(P_ -> -minimumP_helper(k0, kin, s, P_, N_points, err_points, E_quadrature, E_multipole), -tol, P_min, P_max) @@ -92,7 +92,7 @@ end function minimumP_helper(k0, kin, s, P, N_points, err_points, E_quadrature, E_multipole) #helper function for computing error for given P - E_multipole[:] = 0.0 + E_multipole[:] .= 0.0 sp = ScatteringProblem([s],[1],[0.0 0.0],[0.0]) beta_p,inn = solve_particle_scattering(k0, kin, P, sp, PlaneWave(0.0), verbose = false) diff --git a/src/multipole.jl b/src/multipole.jl index 627204f..b0704d8 100644 --- a/src/multipole.jl +++ b/src/multipole.jl @@ -13,62 +13,63 @@ printed if `verbose` is true. """ function solve_particle_scattering(k0, kin, P, sp::ScatteringProblem, u::Einc; get_inner = true, verbose = true) + #TODO: incorporate matrix-lexx rotation rotateMultipole! shapes = sp.shapes; ids = sp.ids; centers = sp.centers; φs = sp.φs Ns = size(sp) #first solve for single scatterer densities - tic() - scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids) - dt1 = toq() - tic() - T = Array{Complex{Float64}}(Ns*(2*P+1),Ns*(2*P+1)) - a = u2α(k0, u, centers, P) - Btemp = Array{Complex{Float64}}(2*P+1,2*P+1) - for ic1 = 1:Ns - rng1 = (ic1-1)*(2*P+1) + (1:2*P+1) - if φs[ic1] == 0.0 || typeof(shapes[ids[ic1]]) == CircleParams - RotScatMat = scatteringMatrices[ids[ic1]] - else - Rot = spdiagm(Complex{Float64}[exp(-1.0im*φs[ic1]*l) for l=-P:P]) #rotation matrix - RotScatMat = Rot*(scatteringMatrices[ids[ic1]]*conj(Rot)) - end - for ic2 = 1:Ns - rng2 = (ic2-1)*(2*P+1) + (1:2*P+1) - if ic1 == ic2 - T[rng1,rng2] = eye(Complex{Float64},2*P+1) + dt1 = @elapsed begin + scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids) + end + dt2 = @elapsed begin + T = Array{Complex{Float64}}(undef, Ns*(2*P+1), Ns*(2*P+1)) + a = u2α(k0, u, centers, P) + Btemp = Array{Complex{Float64}}(undef, 2*P+1, 2*P+1) + for ic1 = 1:Ns + rng1 = (ic1-1)*(2*P+1) .+ (1:2*P+1) + if φs[ic1] == 0.0 || typeof(shapes[ids[ic1]]) == CircleParams + RotScatMat = scatteringMatrices[ids[ic1]] else - M2Lmatrix!(Btemp, k0, P, centers[ic1,:] - centers[ic2,:]) - T[rng1,rng2] = -RotScatMat*Btemp + Rot = sparse(Diagonal(Complex{Float64}[exp(-1.0im*φs[ic1]*l) for l=-P:P])) #rotation matrix + RotScatMat = Rot*(scatteringMatrices[ids[ic1]]*conj(Rot)) + end + for ic2 = 1:Ns + rng2 = (ic2-1)*(2*P+1) .+ (1:2*P+1) + if ic1 == ic2 + T[rng1,rng2] = Matrix{Complex{Float64}}(I, 2*P+1, 2*P+1) + else + M2Lmatrix!(Btemp, k0, P, centers[ic1,:] - centers[ic2,:]) + T[rng1,rng2] = -RotScatMat*Btemp + end end + a[rng1] = RotScatMat*a[rng1] end - a[rng1] = RotScatMat*a[rng1] end - dt2 = toq() - tic() - beta = T\a - dt3 = toq() + dt3 = @elapsed begin + beta = T\a + end #recover full incoming expansion - in sigma_mu terms for parametrized shape, #in multipole expansion for circle if get_inner - tic() #find LU factorization once for each shape - scatteringLU = [lufact(scatteringMatrices[iid]) for iid = 1:length(shapes)] - inner = Array{Vector{Complex{Float64}}}(Ns) - α_c = Array{Complex{Float64}}(2*P+1) - for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) - if typeof(shapes[ids[ic]]) == ShapeParams - if φs[ic] == 0.0 - α_c[:] = scatteringLU[ids[ic]]\beta[rng] + dt4 = @elapsed begin + scatteringLU = [lu(scatteringMatrices[iid]) for iid = 1:length(shapes)] + inner = Array{Vector{Complex{Float64}}}(undef, Ns) + α_c = Array{Complex{Float64}}(undef, 2*P+1) + for ic = 1:Ns + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) + if typeof(shapes[ids[ic]]) == ShapeParams + if φs[ic] == 0.0 + α_c[:] = scatteringLU[ids[ic]]\beta[rng] + else + Rot = sparse(Diagonal(Complex{Float64}[exp(-1.0im*φs[ic]*l) for l=-P:P])) #rotation matrix + α_c[:] = scatteringLU[ids[ic]]\(conj(Rot)*beta[rng]) + end + inner[ic] = innerExpansions[ids[ic]]*α_c else - Rot = spdiagm(Complex{Float64}[exp(-1.0im*φs[ic]*l) for l=-P:P]) #rotation matrix - α_c[:] = scatteringLU[ids[ic]]\(conj(Rot)*beta[rng]) + inner[ic] = innerExpansions[ids[ic]]*beta[rng] end - inner[ic] = innerExpansions[ids[ic]]*α_c - else - inner[ic] = innerExpansions[ids[ic]]*beta[rng] end end - dt4 = toq() end verbose && begin println("Direct solution timing:") @@ -81,18 +82,18 @@ function solve_particle_scattering(k0, kin, P, sp::ScatteringProblem, u::Einc; end function M2Lmatrix!(T, k, P, d) - #builds non-diagonal translation matrix from particle 2 to particle 1 (d=1-2) + #builds off-diagonal translation matrix from particle 2 to particle 1 (d=1-2) #only calculates 2*(2*P+1) values kd = k*sqrt(sum(abs2,d)) - td = atan2(d[2],d[1]) + td = atan(d[2],d[1]) bess = besselh.(0:2*P,1,kd) for ix = 1:2*P #lower diagonals rng = ix+1:1+(2*P+1):(2*P+1)^2-(2*P+1)*ix - T[rng] = exp(-1im*td*ix)*(-1)^(ix)*bess[ix+1] + T[rng] .= exp(-1im*td*ix)*(-1)^(ix)*bess[ix+1] end for ix = 0:2*P #central and upper diagonals rng = ix*(2*P+1)+1:1+(2*P+1):(2*P+1)^2-ix - T[rng] = exp(1im*td*ix)*bess[ix+1] + T[rng] .= exp(1im*td*ix)*bess[ix+1] end end @@ -118,7 +119,7 @@ function scattered_field_multipole!(Esc::Array{Complex{Float64},1}, k0, beta, P, points_moved2 = points[ip,2] - centers[ic,2] rs_moved = hypot(points_moved1, points_moved2) - ts_moved = atan2(points_moved2, points_moved1) + ts_moved = atan(points_moved2, points_moved1) if rs_moved == 0 error("rs_moved == 0, center=$(centers[ic,:]), point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip") end @@ -138,7 +139,7 @@ function scattered_field_multipole_recurrence!(Esc::Array{Complex{Float64},1}, k points_moved2 = points[ip,2] - centers[ic,2] rs_moved = hypot(points_moved1, points_moved2) - ts_moved = atan2(points_moved2, points_moved1) + ts_moved = atan(points_moved2, points_moved1) if rs_moved == 0 error("rs_moved == 0, center=$(centers[ic,:]), point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip") end @@ -156,8 +157,8 @@ end function circleScatteringMatrix(kout, kin, R, P; gamma = false) #non-vectorized, reuses bessel - S = Array{Complex{Float64}}(2*P+1) - gamma && (G = Array{Complex{Float64}}(2*P+1)) + S = Array{Complex{Float64}}(undef, 2*P+1) + gamma && (G = Array{Complex{Float64}}(undef, 2*P+1)) pre_J0 = besselj(-P-1,kout*R) pre_J1 = besselj(-P-1,kin*R) @@ -182,9 +183,9 @@ function circleScatteringMatrix(kout, kin, R, P; gamma = false) pre_J1 = J1 pre_H = H end - ScatMat = spdiagm(S) + ScatMat = sparse(Diagonal(S)) if gamma - InnerMat = spdiagm(G) + InnerMat = sparse(Diagonal(G)) return ScatMat,InnerMat else return ScatMat @@ -197,11 +198,11 @@ function innerFieldCircle(kin, gamma, center::Array{Float64,1}, points::Array{Fl len = size(points,1) points_moved = similar(points) - points_moved[:,1] = points[:,1] - center[1] - points_moved[:,2] = points[:,2] - center[2] + points_moved[:,1] = points[:,1] .- center[1] + points_moved[:,2] = points[:,2] .- center[2] - rs_moved = sqrt.(sum(abs2,points_moved,2)) - ts_moved = atan2.(points_moved[:,2], points_moved[:,1]) + rs_moved = sqrt.(sum(abs2,points_moved, dims=2)) + ts_moved = atan.(points_moved[:,2], points_moved[:,1]) bess = [besselj(p,kin*rs_moved[ii]) for ii=1:len, p=0:P] E = gamma[P + 1]*bess[:,1] @@ -218,7 +219,7 @@ function innerFieldCircle(kin, gamma, center::Array{Float64,1}, points::Array{Fl points_moved = points - center rs_moved = sqrt(sum(abs2,points_moved)) - ts_moved = atan2(points_moved[2], points_moved[1]) + ts_moved = atan(points_moved[2], points_moved[1]) bess = [besselj(p,kin*rs_moved) for p=0:P] E = gamma[P + 1]*bess[1] @@ -229,13 +230,14 @@ function innerFieldCircle(kin, gamma, center::Array{Float64,1}, points::Array{Fl end function particleExpansion(k0, kin, shapes, P, ids) - scatteringMatrices = Array{Any}(0) - innerExpansions = Array{Any}(0) + scatteringMatrices = Array{Any}(undef, 0) + innerExpansions = Array{Any}(undef, 0) for i = 1:length(shapes) #no use in computing matrices if shape doesn't actually show up! push garbage to maintain order + #TODO: don't even try to compute LU for these if all(ids .!= i) - push!(scatteringMatrices, [0.0im 0.0im]) - push!(innerExpansions, [0.0im 0.0im]) + push!(scatteringMatrices, [NaN NaN]) + push!(innerExpansions, [NaN NaN]) continue end if typeof(shapes[i]) == ShapeParams diff --git a/src/optimize_phis.jl b/src/optimize_phis.jl index a6204e3..64d671d 100644 --- a/src/optimize_phis.jl +++ b/src/optimize_phis.jl @@ -48,11 +48,11 @@ end function optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, ids, mFMM, opt, buf) if φs != last_φs - copy!(last_φs, φs) + copyto!(last_φs, φs) #do whatever common calculations and save to shared_var #construct rhs for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) if φs[ic] == 0.0 buf.rhs[rng] = scatteringMatrices[ids[ic]]*α[rng] else @@ -63,19 +63,11 @@ function optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, end end - if opt.method == "pre" - MVP = LinearMap{eltype(buf.rhs)}( - (output_, x_) -> FMM_mainMVP_pre!(output_, x_, - scatteringMatrices, φs, ids, P, mFMM, - buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - elseif opt.method == "pre2" - MVP = LinearMap{eltype(buf.rhs)}( - (output_, x_) -> FMM_mainMVP_pre2!(output_, x_, - scatteringMatrices, φs, ids, P, mFMM, - buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - end + MVP = LinearMap{eltype(buf.rhs)}( + (output_, x_) -> FMM_mainMVP_pre!(output_, x_, + scatteringMatrices, φs, ids, P, mFMM, + buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) fill!(shared_var.β,0.0) shared_var.β,ch = gmres!(shared_var.β, MVP, buf.rhs, @@ -83,7 +75,7 @@ function optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, initially_zero = true) #no restart, preconditioning !ch.isconverged && error("FMM process did not converge") - shared_var.f[:] = H.'*shared_var.β + shared_var.f[:] = transpose(H)*shared_var.β shared_var.f .+= uinc_ end end @@ -99,32 +91,24 @@ function optimize_φ_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, k0, centers,scatteringMatrices, ids, mFMM, opt, buf) - if opt.method == "pre" - MVP = LinearMap{eltype(buf.rhs)}( - (output_, x_) -> FMM_mainMVP_pre!(output_, x_, - scatteringMatrices, φs, ids, P, mFMM, - buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - elseif opt.method == "pre2" - MVP = LinearMap{eltype(buf.rhs)}( - (output_, x_) -> FMM_mainMVP_pre2!(output_, x_, - scatteringMatrices, φs, ids, P, mFMM, - buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - end + MVP = LinearMap{eltype(buf.rhs)}( + (output_, x_) -> FMM_mainMVP_pre!(output_, x_, scatteringMatrices, + φs, ids, P, mFMM, buf.pre_agg, + buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) #time for gradient shared_var.rhs_grad[:] = 0.0 shared_var.∂β[:] = 0.0 D = -1.0im*collect(-P:P) - tempn = Array{Complex{Float64}}(2*P+1) + tempn = Array{Complex{Float64}}(undef, 2*P+1) for n = 1:Ns #compute n-th gradient - rng = (n-1)*(2*P+1) + (1:2*P+1) + rng = (n-1)*(2*P+1) .+ (1:2*P+1) rotateMultipole!(tempn, shared_var.β[rng], -φs[n], P) tempn[:] = scatteringLU[ids[n]]\tempn #LU decomp with pivoting tempn[:] .*= -D v = view(shared_var.rhs_grad, rng) - A_mul_B!(v, scatteringMatrices[ids[n]], tempn) + mul!(v, scatteringMatrices[ids[n]], tempn) rotateMultipole!(v, φs[n], P) v[:] += D.*shared_var.β[rng] @@ -132,7 +116,7 @@ function optimize_φ_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, shared_var.rhs_grad, restart = Ns*(2*P+1), tol = 10*opt.tol, log = true, initially_zero = true) - #warn("using dbdn_tol = 10*opt.tol = $(10*opt.tol)") + #@warn("using dbdn_tol = 10*opt.tol = $(10*opt.tol)") if ch.isconverged == false display("FMM process did not converge for partial derivative $n/$Ns. ") @@ -143,7 +127,7 @@ function optimize_φ_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, end grad_stor[:] = ifelse(minimize, 2, -2)* - real(shared_var.∂β.'*(H*conj(shared_var.f))) + real(transpose(shared_var.∂β)*(H*conj(shared_var.f))) end @@ -172,10 +156,10 @@ function optimize_φ_adj_g!(grad_stor, φs, shared_var, last_φs, α, H, points, #note: this assumes that there are exactly 2P+1 non zero elements in ∂X. If #not, v must be (2P+1)Ns × 1. D = -1.0im*collect(-P:P) - v = Array{Complex{Float64}}(2P+1) #TODO: minimize dynamic alloc + v = Array{Complex{Float64}}(undef, 2P+1) #TODO: minimize dynamic alloc for n = 1:Ns #compute n-th element of gradient - rng = (n-1)*(2*P+1) + (1:2*P+1) + rng = (n-1)*(2*P+1) .+ (1:2*P+1) rotateMultipole!(v, view(shared_var.β,rng), -φs[n], P) v[:] = scatteringLU[ids[n]]\v #LU decomp with pivoting v[:] .*= -D @@ -183,19 +167,19 @@ function optimize_φ_adj_g!(grad_stor, φs, shared_var, last_φs, α, H, points, rotateMultipole!(v, φs[n], P) v[:] += D.*shared_var.β[rng] - grad_stor[n] = ifelse(minimize, -2, 2)*real(view(λadj,rng).'*v) + grad_stor[n] = ifelse(minimize, -2, 2)*real(transpose(view(λadj,rng))*v) #prepare for next one - #TODO: check why this is here - v[:] = 0.0im + v[:] .= 0.0im end end function optimizationHmatrix(points, centers, Ns, P, k0) - points_moved = Array{Float64}(2) - H = Array{Complex{Float64}}(Ns*(2*P+1), size(points,1)) + points_moved = Array{Float64}(undef, 2) + H = Array{Complex{Float64}}(undef, Ns*(2*P+1), size(points,1)) for ic = 1:Ns, i = 1:size(points,1) points_moved[1] = points[i,1] - centers[ic,1] points_moved[2] = points[i,2] - centers[ic,2] - r_angle = atan2(points_moved[2], points_moved[1]) + r_angle = atan(points_moved[2], points_moved[1]) kr = k0*hypot(points_moved[1], points_moved[2]) ind = (ic-1)*(2*P+1) + P + 1 @@ -216,7 +200,7 @@ function prepare_fmm_reusal_φs(k0, kin, P, shapes, centers, ids, fmmopt) P2, Q = FMMtruncation(fmmopt.acc, boxSize, k0) mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri = true) scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids) - scatteringLU = [lufact(scatteringMatrices[iid]) for iid = 1:length(shapes)] + scatteringLU = [lu(scatteringMatrices[iid]) for iid = 1:length(shapes)] buf = FMMbuffer(Ns,P,Q,length(groups)) return mFMM, scatteringMatrices, scatteringLU, buf end diff --git a/src/optimize_phis_pwr.jl b/src/optimize_phis_pwr.jl index 9e94707..8bad70c 100644 --- a/src/optimize_phis_pwr.jl +++ b/src/optimize_phis_pwr.jl @@ -1,11 +1,11 @@ #for now just the main functions, no wrapper function optimize_pwr_φ_common!(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, Ns, fmmbuf, fmmopts, mFMM) if φs != last_φs - copy!(last_φs, φs) + copyto!(last_φs, φs) #do whatever common calculations and save to shared_var for i in eachindex(opb) for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) if φs[ic] == 0.0 fmmbuf.rhs[rng] = scatteringMatrices[i][ids[ic]]*opb[i].α[rng] else @@ -22,7 +22,7 @@ function optimize_pwr_φ_common!(φs, last_φs, opb, power_buffer, ids, scatteri fmmbuf.pre_agg, fmmbuf.trans), Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - opb[i].β[:] = 0 + opb[i].β[:] .= 0 opb[i].β,ch = gmres!(opb[i].β, MVP, fmmbuf.rhs, restart = Ns*(2*P+1), tol = fmmopts.tol, log = true, initially_zero = true) #no restart @@ -58,7 +58,7 @@ function optimize_pwr_φ_g!(grad_stor, φs, last_φs, opb, power_buffer, ids, sc Ns*(2*P+1), Ns*(2*P+1), ismutating = true) #solve adjoint problem - opb[i].λadj[:] = 0 + opb[i].λadj[:] .= 0 opb[i].λadj, ch = gmres!(opb[i].λadj, MVP, opb[i].rhs_grad, restart = Ns*(2*P+1), tol = fmmopts.tol, log=true, initially_zero = true) #no restart @@ -66,25 +66,19 @@ function optimize_pwr_φ_g!(grad_stor, φs, last_φs, opb, power_buffer, ids, sc normalized residual: $(norm(MVP*opb[i].λadj - opb[i].rhs_grad)/norm(opb[i].rhs_grad))""") D = -1.0im*collect(-P:P) - v = Array{Complex{Float64}}(2*P+1) #TODO: minimize dynamic alloc + v = Array{Complex{Float64}}(undef, 2*P+1) #TODO: minimize dynamic alloc + v2 = Array{Complex{Float64}}(undef, 2*P+1) for n = 1:Ns - rng = (n-1)*(2*P+1) + (1:2*P+1) - # rotateMultipole!(v, view(opb[i].β,rng), -φs[n], P) - # v[:] = scatteringLU[i][ids[n]]\v #LU decomp with pivoting - # v[:] .*= -D - # v[:] = scatteringMatrices[i][ids[n]]*v - # rotateMultipole!(v, φs[n], P) - # v[:] += D.*opb[i].β[rng] - # grad_stor[n] += (-2)*real(view(opb[i].λadj,rng).'*v) - # #prepare for next one - #TODO: check why this is here - # v[:] = 0.0im - - #dumb way - Φ = spdiagm(exp.(-1.0im*(-P:P)*φs[n])) - temppp = Φ*(scatteringLU[i][ids[n]]\(conj(Φ)*opb[i].β[rng])) - temppp2 = spdiagm(D)*(Φ*(scatteringMatrices[i][ids[n]]*(conj(Φ)*temppp))) - - Φ*(scatteringMatrices[i][ids[n]]*(conj(Φ)*(spdiagm(D)*temppp))) - grad_stor[n] += (-2)*real(opb[i].λadj[rng].'*temppp2) + rng = (n-1)*(2*P+1) .+ (1:2*P+1) + rotateMultipole!(v, view(opb[i].β,rng), -φs[n], P) + ldiv!(scatteringLU[i][ids[n]], v) + v2[:] = scatteringMatrices[i][ids[n]]*v + v2 .*= D + v .*= D + v[:] = scatteringMatrices[i][ids[n]]*v + v2 .-= v + rotateMultipole!(v2, φs[n], P) + grad_stor[n] += (-2)*real(transpose(opb[i].λadj[rng])*v2) end end end diff --git a/src/optimize_rs.jl b/src/optimize_rs.jl index d82da93..78f5bcc 100644 --- a/src/optimize_rs.jl +++ b/src/optimize_rs.jl @@ -23,16 +23,16 @@ function optimize_radius(rs0, r_min, r_max, points, ids, P, ui, k0, kin, Ns = size(centers,1) J = length(rs0) - assert(maximum(ids) <= J) + @assert maximum(ids) <= J if length(r_min) == 1 r_min = r_min*ones(Float64,J) else - assert(J == length(r_min)) + @assert J == length(r_min) end if length(r_max) == 1 r_max = r_max*ones(Float64,J) else - assert(J == length(r_max)) + @assert J == length(r_max) end verify_min_distance(CircleParams.(r_max), centers, ids, points) || error("Particles are too close or r_max are too large.") @@ -43,8 +43,8 @@ function optimize_radius(rs0, r_min, r_max, points, ids, P, ui, k0, kin, mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri = true) #allocate derivative - scatteringMatrices = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] - dS_S = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] + scatteringMatrices = [sparse(one(Complex{Float64})I, 2*P+1, 2*P+1) for ic = 1:J] + dS_S = [sparse(one(Complex{Float64})I, 2*P+1, 2*P+1) for ic = 1:J] #stuff that is done once H = optimizationHmatrix(points, centers, Ns, P, k0) @@ -56,7 +56,7 @@ function optimize_radius(rs0, r_min, r_max, points, ids, P, ui, k0, kin, buf = FMMbuffer(Ns,P,Q,length(groups)) shared_var = OptimBuffer(Ns,P,size(points,1),J) initial_rs = copy(rs0) - last_rs = similar(initial_rs); last_rs[1] = NaN; assert(last_rs != initial_rs) #initial_rs, last_rs must be different before first iteration + last_rs = similar(initial_rs); last_rs[1] = NaN; @assert last_rs != initial_rs #initial_rs, last_rs must be different before first iteration if adjoint gopt! = (grad_stor, rs) -> optimize_radius_adj_g!(grad_stor, rs, last_rs, @@ -86,31 +86,25 @@ end function optimize_radius_common!(rs, last_rs, shared_var, φs, α, H, points, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf) if (rs != last_rs) - copy!(last_rs, rs) + copyto!(last_rs, rs) #do whatever common calculations and save to shared_var #construct rhs for id in unique(ids) try updateCircleScatteringDerivative!(scatteringMatrices[id], dS_S[id], k0, kin, rs[id], P) catch - warn("Could not calculate derivatives for id=$id,k0=$k0,kin=$kin,R=$(rs[id])") + @warn("Could not calculate derivatives for id=$id,k0=$k0,kin=$kin,R=$(rs[id])") error() end end for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) buf.rhs[rng] = scatteringMatrices[ids[ic]]*α[rng] end - if opt.method == "pre" - MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre!(output_, - x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - elseif opt.method == "pre2" - MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre2!(output_, - x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - end + MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre!(output_, + x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) fill!(shared_var.β,0.0) #no restart @@ -120,7 +114,7 @@ function optimize_radius_common!(rs, last_rs, shared_var, φs, α, H, points, P, if !ch.isconverged error("FMM process did not converge") end - shared_var.f[:] = H.'*shared_var.β + shared_var.f[:] = transpose(H)*shared_var.β shared_var.f .+= uinc_ #incident end end @@ -134,15 +128,9 @@ end function optimize_radius_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, points, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf, minimize) optimize_radius_common!(rs, last_rs, shared_var, φs, α, H, points, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf) - if opt.method == "pre" - MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre!(output_, - x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - elseif opt.method == "pre2" - MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre2!(output_, - x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - end + MVP = LinearMap{eltype(buf.rhs)}((output_, x_) -> FMM_mainMVP_pre!(output_, + x_, scatteringMatrices, φs, ids, P, mFMM, buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) #time for gradient shared_var.∂β[:] = 0.0 @@ -151,7 +139,7 @@ function optimize_radius_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, poin #as more than one beta is affected. #TODO:expand to phi optimization?, and clean+speed this up. for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) if ids[ic] == n shared_var.rhs_grad[rng] = dS_S[n]*shared_var.β[rng] else @@ -173,7 +161,7 @@ function optimize_radius_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, poin end end - grad_stor[:] = ifelse(minimize,2,-2)*real(shared_var.∂β.'*(H*conj(shared_var.f))) + grad_stor[:] = ifelse(minimize,2,-2)*real(transpose(shared_var.∂β)*(H*conj(shared_var.f))) end function optimize_radius_adj_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, points, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf, minimize) @@ -200,16 +188,17 @@ function optimize_radius_adj_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, for n = 1:length(rs) #compute n-th gradient - here we must pay the price for symmetry #as more than one beta is affected. Overwrites rhs_grad. + #TODO: only sum over relevant parts for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) if ids[ic] == n shared_var.rhs_grad[rng] = dS_S[n]*shared_var.β[rng] else - shared_var.rhs_grad[rng] = 0.0 + shared_var.rhs_grad[rng] .= 0.0 end end - grad_stor[n] = ifelse(minimize,-2,2)*real(λadj.'*shared_var.rhs_grad) + grad_stor[n] = ifelse(minimize,-2,2)*real(transpose(λadj)*shared_var.rhs_grad) end end @@ -259,8 +248,8 @@ end # mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri = true) # # #allocate derivative -# scatteringMatrices = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] -# dS_S = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] +# scatteringMatrices = [sparse(one(Complex{Float64})I, 2*P+1, 2*P+1) for ic = 1:J] +# dS_S = [sparse(one(Complex{Float64})I, 2*P+1, 2*P+1) for ic = 1:J] # # #stuff that is done once # H = optimizationHmatrix(points, centers, Ns, P, k0) @@ -273,7 +262,7 @@ end # shared_var = OptimBuffer(Ns,P,size(points,1),J) # last_rs = similar(rs); all(last_rs .== rs) && (last_rs[1] += 1) # -# grad_stor = Array{Float64}(J) +# grad_stor = Array{Float64}(undef, J) # optimize_radius_adj_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, points, # P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, # fmmopts, buf, minimize) diff --git a/src/optimize_rs_mf.jl b/src/optimize_rs_mf.jl index a2a9a1f..dc66e27 100644 --- a/src/optimize_rs_mf.jl +++ b/src/optimize_rs_mf.jl @@ -6,18 +6,18 @@ function optimize_rs_mf(rs0, r_min, r_max, points, point_flags, ids, P, ui, k0, Ns = size(centers,1) φs = zeros(Ns) - Nk = length(k0); assert(Nk == length(kin)) + Nk = length(k0); @assert Nk == length(kin) J = length(rs0) - assert(maximum(ids) <= J) + @assert maximum(ids) <= J if length(r_min) == 1 r_min = r_min*ones(Float64,J) else - assert(J == length(r_min)) + @assert J == length(r_min) end if length(r_max) == 1 r_max = r_max*ones(Float64,J) else - assert(J == length(r_max)) + @assert J == length(r_max) end verify_min_distance(CircleParams.(r_max), centers, ids, points) || error("Particles are too close or r_max are too large.") @@ -25,8 +25,8 @@ function optimize_rs_mf(rs0, r_min, r_max, points, point_flags, ids, P, ui, k0, groups,boxSize = divideSpace(centers, fmmopts) P2,Q = FMMtruncation(fmmopts.acc, boxSize, maximum(k0)) mFMM = [FMMbuildMatrices(k0[ik], P, P2, Q, groups, centers, boxSize) for ik = 1:Nk] - scatteringMatrices = [[speye(Complex{Float64}, 2*P+1) for i = 1:J] for ik = 1:Nk] - dS_S = [[speye(Complex{Float64}, 2*P+1) for i = 1:J] for ik = 1:Nk] + scatteringMatrices = [[sparse(one(Complex{Float64})I, 2*P+1, 2*P+1) for i = 1:J] for ik = 1:Nk] + dS_S = [[sparse(one(Complex{Float64})I, 2*P+1, 2*P+1) for i = 1:J] for ik = 1:Nk] #calculate nd expand incident field α = [u2α(k0[ik], ui, centers, P) for ik = 1:Nk] @@ -37,7 +37,7 @@ function optimize_rs_mf(rs0, r_min, r_max, points, point_flags, ids, P, ui, k0, buf = FMMbuffer(Ns,P,Q,length(groups)) shared_var = [OptimBuffer(Ns, P, size(points,1), J) for ik = 1:Nk] initial_rs = copy(rs0) - last_rs = similar(initial_rs); last_rs[1] = NaN; assert(last_rs != initial_rs) #initial_rs, last_rs must be different before first iteration + last_rs = similar(initial_rs); last_rs[1] = NaN; @assert last_rs != initial_rs #initial_rs, last_rs must be different before first iteration df = OnceDifferentiable( rs -> optimize_rs_mf_f(rs, last_rs, shared_var, φs, α, H, points, point_flags, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, @@ -54,7 +54,7 @@ function optimize_rs_mf_common!(rs, last_rs, shared_var, φs, α, H, points, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, fmmopts, buf) if rs != last_rs - copy!(last_rs, rs) + copyto!(last_rs, rs) #do whatever common calculations and save to shared_var #construct rhs for ik in eachindex(k0) @@ -63,7 +63,7 @@ function optimize_rs_mf_common!(rs, last_rs, shared_var, φs, α, dS_S[ik][id], k0[ik], kin[ik], rs[id], P) end for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) buf.rhs[rng] = scatteringMatrices[ik][ids[ic]]*α[ik][rng] end @@ -80,7 +80,7 @@ function optimize_rs_mf_common!(rs, last_rs, shared_var, φs, α, !ch.isconverged && error("""FMM process did not converge, normalized residual: $(norm(MVP*shared_var[ik].β - buf.rhs)/norm(buf.rhs))""") - shared_var[ik].f[:] = H[ik].'*shared_var[ik].β + shared_var[ik].f[:] = transpose(H[ik])*shared_var[ik].β shared_var[ik].f[:] .+= uinc_[ik] end end @@ -121,14 +121,14 @@ function optimize_rs_mf_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, point #compute n-th gradient - here we must pay the price for symmetry #as more than one beta is affected. Overwrites rhs_grad. for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) if ids[ic] == n shared_var[ik].rhs_grad[rng] = dS_S[ik][n]*shared_var[ik].β[rng] else shared_var[ik].rhs_grad[rng] = 0.0 end end - grad_stor[n] += -2*real(λadj.'*shared_var[ik].rhs_grad) + grad_stor[n] += -2*real(transpose(λadj)*shared_var[ik].rhs_grad) end end end diff --git a/src/optimize_rs_pwr.jl b/src/optimize_rs_pwr.jl index c0ff044..b8a3c3c 100644 --- a/src/optimize_rs_pwr.jl +++ b/src/optimize_rs_pwr.jl @@ -18,19 +18,19 @@ mutable struct PowerBuffer #one for each *power calculation* (can be > 1 for eac l::Float64 #arc length PowerBuffer(Ns::Integer, P::Integer, nhat, iₚ) = new(iₚ, - Array{Complex{Float64}}(size(nhat, 1)), - Array{Complex{Float64}}(size(nhat, 1)), - Array{Complex{Float64}}(size(nhat, 1)), - Array{Complex{Float64}}(size(nhat, 1)), - Array{Complex{Float64}}(size(nhat, 1)), - Array{Complex{Float64}}(size(nhat, 1)), + Array{Complex{Float64}}(undef, size(nhat, 1)), + Array{Complex{Float64}}(undef, size(nhat, 1)), + Array{Complex{Float64}}(undef, size(nhat, 1)), + Array{Complex{Float64}}(undef, size(nhat, 1)), + Array{Complex{Float64}}(undef, size(nhat, 1)), + Array{Complex{Float64}}(undef, size(nhat, 1)), 0.0, - [Array{Complex{Float64}}(Ns*(2P+1)) for i=1:size(nhat, 1)], - [Array{Complex{Float64}}(Ns*(2P+1)) for i=1:size(nhat, 1)], - [Array{Complex{Float64}}(Ns*(2P+1)) for i=1:size(nhat, 1)], + [Array{Complex{Float64}}(undef, Ns*(2P+1)) for i=1:size(nhat, 1)], + [Array{Complex{Float64}}(undef, Ns*(2P+1)) for i=1:size(nhat, 1)], + [Array{Complex{Float64}}(undef, Ns*(2P+1)) for i=1:size(nhat, 1)], nhat, - Array{Complex{Float64},2}(size(nhat)), - Array{Complex{Float64}}(Ns*(2P+1)), + Array{Complex{Float64},2}(undef, size(nhat)), + Array{Complex{Float64}}(undef, Ns*(2P+1)), 0.0) end @@ -38,15 +38,15 @@ end function optMatrixPwr(points, centers, Ns, P, k0, ui, iₚ, nhat, l) sv = PowerBuffer(Ns, P, nhat, iₚ) cf = 1/(1im*k0*eta0) - pt = Array{Float64}(2) + pt = Array{Float64}(undef, 2) for ip = 1:size(points,1) for ic = 1:Ns pt[1] = points[ip,1] - centers[ic,1] pt[2] = points[ip,2] - centers[ic,2] - θ = atan2(pt[2], pt[1]) + θ = atan(pt[2], pt[1]) R = hypot(pt[1], pt[2]) - ind = (ic-1)*(2*P+1) + P + 1 + ind = (ic-1)*(2*P+1) .+ P + 1 Hₚ₋₂ = besselh(-1, 1, k0*R) Hₚ₋₁ = besselh(0, 1, k0*R) @@ -88,13 +88,15 @@ end function OptimProblemBuffer(k0::Float64, kin::Float64, centers::Array{Float64,2}, ui::Einc, P::Integer) Ns = size(centers, 1) α = ParticleScattering.u2α(k0, ui, centers, P) - OptimProblemBuffer(k0, kin, α, ui, Array{Complex{Float64}}(Ns*(2*P+1)), - Array{Complex{Float64}}(Ns*(2*P+1)), Array{Complex{Float64}}(Ns*(2*P+1))) + OptimProblemBuffer(k0, kin, α, ui, + Array{Complex{Float64}}(undef, Ns*(2*P+1)), + Array{Complex{Float64}}(undef, Ns*(2*P+1)), + Array{Complex{Float64}}(undef, Ns*(2*P+1))) end function optimize_pwr_rs_common!(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM) if rs != last_rs - copy!(last_rs, rs) + copyto!(last_rs, rs) #do whatever common calculations and save to shared_var for i in eachindex(opb) #TODO: if multiple problems have same k0, these calcs are duplicated @@ -102,7 +104,7 @@ function optimize_pwr_rs_common!(rs, last_rs, opb, power_buffer, ids, scattering ParticleScattering.updateCircleScatteringDerivative!(scatteringMatrices[i][id], dS_S[i][id], opb[i].k0, opb[i].kin, rs[id], P) end for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) fmmbuf.rhs[rng] = scatteringMatrices[i][ids[ic]]*opb[i].α[rng] end @@ -112,7 +114,7 @@ function optimize_pwr_rs_common!(rs, last_rs, opb, power_buffer, ids, scattering fmmbuf.pre_agg, fmmbuf.trans), Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - opb[i].β[:] = 0 + opb[i].β[:] .= 0 opb[i].β,ch = gmres!(opb[i].β, MVP, fmmbuf.rhs, restart = Ns*(2*P+1), tol = fmmopts.tol, log = true, initially_zero = true) #no restart @@ -148,7 +150,7 @@ function optimize_pwr_rs_g!(grad_stor, rs, last_rs, opb, power_buffer, ids, scat Ns*(2*P+1), Ns*(2*P+1), ismutating = true) #solve adjoint problem - opb[i].λadj[:] = 0 + opb[i].λadj[:] .= 0 opb[i].λadj, ch = gmres!(opb[i].λadj, MVP, opb[i].rhs_grad, restart = Ns*(2*P+1), tol = fmmopts.tol, log = true, initially_zero = true) #no restart @@ -159,14 +161,14 @@ function optimize_pwr_rs_g!(grad_stor, rs, last_rs, opb, power_buffer, ids, scat #rhs_grad is still usually sparse - utilize this to reduce complexity #here O(N^2) -> O(N) for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) + rng = (ic-1)*(2*P+1) .+ (1:2*P+1) if ids[ic] == n opb[i].rhs_grad[rng] = dS_S[i][n]*opb[i].β[rng] else - opb[i].rhs_grad[rng] = 0.0 + opb[i].rhs_grad[rng] .= 0.0 end end - grad_stor[n] += -2*real(opb[i].λadj.'*opb[i].rhs_grad) + grad_stor[n] += -2*real(transpose(opb[i].λadj)*opb[i].rhs_grad) end end end @@ -174,12 +176,12 @@ end function calc_multi_pwr!(power_buffer, opb) for sv in power_buffer len = length(sv.HEz) - Sn = Array{Complex{Float64}}(len) # S⋅n + Sn = Array{Complex{Float64}}(undef, len) # S⋅n βi = opb[sv.iₚ].β for ip = 1:len - sv.Ez[ip] = sv.HEz[ip].'*βi + sv.Ez_inc[ip] - sv.Hx[ip] = sv.HHx[ip].'*βi + sv.Hx_inc[ip] - sv.Hy[ip] = sv.HHy[ip].'*βi + sv.Hy_inc[ip] + sv.Ez[ip] = transpose(sv.HEz[ip])*βi + sv.Ez_inc[ip] + sv.Hx[ip] = transpose(sv.HHx[ip])*βi + sv.Hx_inc[ip] + sv.Hy[ip] = transpose(sv.HHy[ip])*βi + sv.Hy_inc[ip] Sn[ip] = -0.5*real(sv.Ez[ip]*conj(sv.Hy[ip]))*sv.nhat[ip,1] Sn[ip] += 0.5*real(sv.Ez[ip]*conj(sv.Hx[ip]))*sv.nhat[ip,2] end @@ -192,7 +194,7 @@ function dPdβ_pwr!(sv::PowerBuffer) fill!(sv.∂pow, 0.0) #this is a sum of real and imaginary parts, hence the additional 1/2 factor for ip = 1:len - cf = ((ip == 1 || ip == len) ? 0.5 : 1.0)*sv.l/(len-1) #trapezoidal rule constant + cf = ifelse(ip == 1 || ip == len, 0.5, 1.0)*sv.l/(len-1) #trapezoidal rule constant sv.∂pow .+= (-0.25*cf*sv.nhat[ip,1])*(conj(sv.Hy[ip])*sv.HEz[ip] + conj(sv.Ez[ip])*sv.HHy[ip]) sv.∂pow .+= (0.25*cf*sv.nhat[ip,2])*(conj(sv.Hx[ip])*sv.HEz[ip] + diff --git a/src/poynting.jl b/src/poynting.jl index f2d9475..17057ab 100644 --- a/src/poynting.jl +++ b/src/poynting.jl @@ -3,18 +3,18 @@ function poynting_vector(k0, beta, centers, points, Ez_inc, Hx_inc, Hy_inc) Ns = size(centers,1) P = div(div(length(beta),Ns)-1,2) len = size(points,1) - pyntg = Array{Float64}(len, 2) + pyntg = Array{Float64}(undef, len, 2) Ez = copy(Ez_inc) for ip in 1:len ∂Ez∂x = 0.0im ∂Ez∂y = 0.0im for ic in 1:Ns - ind = (ic-1)*(2*P+1) + P + 1 + ind = (ic-1)*(2*P+1) .+ P + 1 pt1 = points[ip,1] - centers[ic,1] pt2 = points[ip,2] - centers[ic,2] R = hypot(pt1, pt2) - θ = atan2(pt2, pt1) + θ = atan(pt2, pt1) R == 0 && error("R == 0, center=$(centers[ic,:]), point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip") @@ -51,7 +51,7 @@ end function poynting_vector(k0, points, Ez_inc, Hx_inc, Hy_inc) #assumes points lie outside all scattering disks len = size(points,1) - pyntg = Array{Float64}(len, 2) + pyntg = Array{Float64}(undef, len, 2) for ip in 1:len pyntg[ip,1] = (-0.5)*real(Ez_inc[ip]*conj(Hy_inc[ip])) pyntg[ip,2] = 0.5*real(Ez_inc[ip]*conj(Hx_inc[ip])) @@ -61,7 +61,7 @@ end function calc_power(k0::Array, kin, P, sp, points, nhat, ui) #this assumes points are equidistant. correct result only after multiplying by arc length - power = Array{Float64}(length(k0)) + power = Array{Float64}(undef, length(k0)) for i in eachindex(k0) Ez_inc = uinc(k0[i], points, ui) Hx_inc = hxinc(k0[i], points, ui) @@ -102,10 +102,10 @@ end function rect_border(b, Nx, Ny) #clockwise from top - points1 = [linspace(b[1], b[2], Nx) b[4]*ones(Nx)] - points2 = [b[2]*ones(Ny) linspace(b[4], b[3], Ny)] - points3 = [linspace(b[2], b[1], Nx) b[3]*ones(Nx)] - points4 = [b[1]*ones(Ny) linspace(b[3], b[4], Ny)] + points1 = [range(b[1], stop=b[2], length=Nx) b[4]*ones(Nx)] + points2 = [b[2]*ones(Ny) range(b[4], stop=b[3], length=Ny)] + points3 = [range(b[2], stop=b[1], length=Nx) b[3]*ones(Nx)] + points4 = [b[1]*ones(Ny) range(b[3], stop=b[4], length=Ny)] n = [[zeros(Nx) ones(Nx)], [ones(Ny) zeros(Ny)], [zeros(Nx) -ones(Nx)], [-ones(Ny) zeros(Ny)]] [points1, points2, points3, points4], n diff --git a/src/scattering.jl b/src/scattering.jl index 2ce0796..06e4ab7 100644 --- a/src/scattering.jl +++ b/src/scattering.jl @@ -5,71 +5,38 @@ Given a shape `s` with `2N` discretization nodes, outer and inner wavenumbers `kout`,`kin`, and the cylindrical harmonics parameter `P`, returns the potential densities `sigma_mu`. Each column contains the response to a different harmonic, where the first `2N` entries contain the single-layer potential density -(``\sigma``), and the lower entries contain the double-layer density (``\mu``). +(``\\sigma``), and the lower entries contain the double-layer density (``\\mu``). """ function get_potential(kout, kin, P, t, ft, dft) N = length(t) #N here is 2N elesewhere. A = SDNTpotentialsdiff(kout, kin, t, ft, dft) - LU = lufact(A) + LU = lu(A) - sigma_mu = Array{Complex{Float64}}(2*N, 2*P+1) + sigma_mu = Array{Complex{Float64}}(undef, 2*N, 2*P+1) #assuming the wave is sampled on the shape - nz = sqrt.(sum(abs2,ft,2)) - θ = atan2.(ft[:,2], ft[:,1]) - ndz = sqrt.(sum(abs2,dft,2)) + nz = sqrt.(sum(abs2, ft, dims=2)) + θ = atan.(ft[:,2], ft[:,1]) + ndz = sqrt.(sum(abs2, dft, dims=2)) nzndz = nz.*ndz wro = dft[:,2].*ft[:,1] - dft[:,1].*ft[:,2] zz = dft[:,1].*ft[:,1] + dft[:,2].*ft[:,2] bessp = besselj.(-P-1, kout*nz) bess = similar(bessp) - du = Array{Complex{Float64}}(length(bessp)) - rhs = Array{Complex{Float64}}(2*length(bessp)) + du = Array{Complex{Float64}}(undef, length(bessp)) + rhs = Array{Complex{Float64}}(undef, 2*length(bessp)) for p = -P:P bess[:] = besselj.(p, kout*nz) du[:] = kout*bessp.*wro - (p*bess./nz).*(wro + 1im*zz) rhs[:] = -[bess.*exp.(1.0im*p*θ); (du./nzndz).*exp.(1.0im*p*θ)] sigma_mu[:,p + P + 1] = LU\rhs - copy!(bessp, bess) + copyto!(bessp, bess) end return sigma_mu end -# -# function get_potential_deprecated(kout, kin, P, t, ft, dft) -# #problematic assumption here that t=θ in the Jacobi-Anger expansion (for rhs). -# #this means that we assume all shapes are given in polar form (r(θ)cosθ,r(θ)sinθ)), -# #which is not the case for the simple ellipse (r1cosθ,r2sinθ). -# N = length(t) #N here is 2N elesewhere. -# -# A = SDNTpotentialsdiff(kout, kin, t, ft, dft) -# LU = lufact(A) -# -# sigma_mu = Array{Complex{Float64}}(2*N, 2*P+1) -# -# #assuming the wave is sampled on the shape -# nz = sqrt.(sum(abs2,ft,2)) -# ndz = sqrt.(sum(abs2,dft,2)) -# nzndz = nz.*ndz -# wro = dft[:,2].*ft[:,1] - dft[:,1].*ft[:,2] -# zz = dft[:,1].*ft[:,1] + dft[:,2].*ft[:,2] -# -# bessp = besselj.(-P-1, kout*nz) -# bess = similar(bessp) -# du = Array{Complex{Float64}}(length(bessp)) -# rhs = Array{Complex{Float64}}(2*length(bessp)) -# for p = -P:P -# bess[:] = besselj.(p, kout*nz) -# du[:] = kout*bessp.*wro - (p*bess./nz).*(wro + 1im*zz) -# rhs[:] = -[bess.*exp.(1.0im*p*t); -# (du./nzndz).*exp.(1.0im*p*t)] -# sigma_mu[:,p + P + 1] = LU\rhs -# copy!(bessp, bess) -# end -# return sigma_mu -# end """ get_potential(kout, kin, P, t, ft, dft) -> sigma_mu @@ -84,14 +51,15 @@ function SDNTpotentialsdiff(k1, k2, t, ft, dft) (N = div(length(t),2)) : (error("length(t) must be even")) (Rvec, kmlogvec) = KM_weights(N) - A = Array{Complex{Float64}}(4*N, 4*N) + A = Array{Complex{Float64}}(undef, 4*N, 4*N) - ndft = sqrt.(vec(sum(abs2,dft,2))) - rij = Array{Float64}(2) + ndft = sqrt.(vec(sum(abs2,dft,dims=2))) + rij = Array{Float64}(undef, 2) for i=1:2*N, j=1:i if i == j T1 = -(k1^2 - k2^2) - T2 = k1^2*(π*1im - 2γ + 1 - 2*log(k1*ndft[i]/2)) - k2^2*(π*1im - 2γ + 1 - 2*log(k2*ndft[i]/2)) + T2 = k1^2*(π*1im - 2MathConstants.γ + 1 - 2*log(k1*ndft[i]/2)) - + k2^2*(π*1im - 2MathConstants.γ + 1 - 2*log(k2*ndft[i]/2)) A[i,j] = (-ndft[i]/2/N)*log(k1/k2) #dS (dM1=0) A[i,j+2*N] = 1 #dD (dL=0) @@ -156,7 +124,7 @@ function KM_weights(N) #Input: N (integer>=1) #Output: R,K (float vectors of length 2N) arg1 = Float64[cos(m*j*π/N)/m for m=1:N-1, j=0:2*N-1] - R = vec((-2π/N)*sum(arg1,1)) - (π/N^2)*Float64[cos(j*π) for j=0:2*N-1] + R = vec((-2π/N)*sum(arg1,dims=1)) - (π/N^2)*Float64[cos(j*π) for j=0:2*N-1] K = Float64[2*log(2*sin(0.5π*j/N)) for j = 0:2*N-1] return (R,K) @@ -168,16 +136,16 @@ end Given a shape `s` with `2N` discretization nodes, outer and inner wavenumbers `kout`,`kin`, and an incident plane-wave angle, returns the potential densities vector `sigma_mu`. The first `2N` entries contain the single-layer -potential density (``\sigma``), and the lower entries contain the double-layer -density (``\mu``). +potential density (``\\sigma``), and the lower entries contain the double-layer +density (``\\mu``). """ function get_potentialPW(kout, kin, s, θ_i) N = length(s.t) #N here is different... A = SDNTpotentialsdiff(kout, kin, s.t, s.ft, s.dft) - LU = lufact(A) + LU = lu(A) - ndft = sqrt.(sum(abs2,s.dft,2)) + ndft = sqrt.(sum(abs2, s.dft, dims=2)) ui = exp.(1.0im*kout*(cos(θ_i)*s.ft[:,1] + sin(θ_i)*s.ft[:,2])) rhs = -[ui; (1.0im*kout*ui).*((cos(θ_i)*s.dft[:,2] - sin(θ_i)*s.dft[:,1])./ndft)] @@ -207,13 +175,13 @@ function scatteredfield(sigma_mu, k, t, ft, dft, p) #in space with wavenumber k at points p *off* the boundary. For field on the boundary, #SDpotentials function must be used. if size(p,2) == 1 #single point, rotate it - p = p.' + p = transpose(p) end N = length(t) M = size(p,1) r = zeros(Float64,2) #loop is faster here: - SDout = Array{Complex{Float64}}(M, 2*N) + SDout = Array{Complex{Float64}}(undef, M, 2*N) for j = 1:N ndft = hypot(dft[j,1],dft[j,2]) for i = 1:M @@ -221,7 +189,7 @@ function scatteredfield(sigma_mu, k, t, ft, dft, p) nr = hypot(r[1],r[2]) if nr < eps() #TODO: use SDNTpotentialsdiff here - warn("Encountered singularity in scatteredfield.") + @warn("Encountered singularity in scatteredfield.") SDout[i,j] = 0 SDout[i,j+N] = 0 continue @@ -237,10 +205,10 @@ end function shapeMultipoleExpansion(k, t, ft, dft, P) #unlike others (so far), this does *not* assume t_j=pi*j/N N = div(length(t),2) - nz = vec(sqrt.(sum(abs2,ft,2))) - θ = atan2.(ft[:,2], ft[:,1]) - ndz = vec(sqrt.(sum(abs2,dft,2))) - AB = Array{Complex{Float64}}(2*P + 1, 4*N) + nz = vec(sqrt.(sum(abs2, ft, dims=2))) + θ = atan.(ft[:,2], ft[:,1]) + ndz = vec(sqrt.(sum(abs2, dft, dims=2))) + AB = Array{Complex{Float64}}(undef, 2*P + 1, 4*N) bessp = besselj.(-P-1,k*nz) bess = similar(bessp) for l = -P:0 @@ -256,50 +224,25 @@ function shapeMultipoleExpansion(k, t, ft, dft, P) AB[l+P+1,j+2*N] = 0.25im*(π/N)*(exp(-1.0im*l*θ[j])/nz[j])*(b1 + b2) l != 0 && (AB[-l+P+1,j+2*N] = 0.25im*((-1.0)^l*π/N)*(exp(1.0im*l*θ[j])/nz[j])*(-b1_ + b2)) end - copy!(bessp,bess) + copyto!(bessp,bess) end return AB end -# -# function shapeMultipoleExpansion_deprecated(k, t, ft, dft, P) -# #unlike others (so far), this does *not* assume t_j=pi*j/N -# N = div(length(t),2) -# nz = vec(sqrt.(sum(abs2,ft,2))) -# ndz = vec(sqrt.(sum(abs2,dft,2))) -# AB = Array{Complex{Float64}}(2*P + 1, 4*N) -# bessp = besselj.(-P-1,k*nz) -# bess = similar(bessp) -# for l = -P:0 -# bess[:] = besselj.(l,k*nz) -# for j = 1:2*N -# AB[l+P+1,j] = 0.25im*(π/N)*bess[j]*exp(-1.0im*l*t[j])*ndz[j] -# l!=0 && (AB[-l+P+1,j] = 0.25im*((-1.0)^l*π/N)*bess[j]*exp(1.0im*l*t[j])*ndz[j]) -# wro = ft[j,1]*dft[j,2] - ft[j,2]*dft[j,1] -# zdz = -1.0im*(ft[j,1]*dft[j,1] + ft[j,2]*dft[j,2]) -# b1 = (-l*bess[j]/nz[j])*(zdz + wro) -# b1_ = (-l*bess[j]/nz[j])*(zdz - wro) -# b2 = k*bessp[j]*wro -# AB[l+P+1,j+2*N] = 0.25im*(π/N)*(exp(-1.0im*l*t[j])/nz[j])*(b1 + b2) -# l!=0 && (AB[-l+P+1,j+2*N] = 0.25im*((-1.0)^l*π/N)*(exp(1.0im*l*t[j])/nz[j])*(-b1_ + b2)) -# end -# copy!(bessp,bess) -# end -# return AB -# end function solvePotential_forError(kin, kout, shape, ls_pos, ls_amp, θ_i) #plane wave outside, line sources inside N = length(shape.t) #N here is different... A = SDNTpotentialsdiff(kout, kin, shape.t, shape.ft, shape.dft) - LU = lufact(A) + LU = lu(A) - ndft = sqrt.(sum(abs2,shape.dft,2)) + ndft = sqrt.(sum(abs2, shape.dft, dims=2)) - r = sqrt.((shape.ft[:,1] - ls_pos[1,1]).^2 + (shape.ft[:,2] - ls_pos[1,2]).^2) + r = sqrt.((shape.ft[:,1] .- ls_pos[1,1]).^2 + + (shape.ft[:,2] .- ls_pos[1,2]).^2) - uls = -ls_amp[1]*0.25im*besselh.(0,kout*r) - duls = ls_amp[1]*0.25im*(kout*besselh.(1,kout*r)).*((shape.ft[:,1]-ls_pos[1,1]).*shape.dft[:,2]-(shape.ft[:,2]-ls_pos[1,2]).*shape.dft[:,1])./r./ndft + uls = (-ls_amp[1]*0.25im)*besselh.(0,kout*r) + duls = (ls_amp[1]*0.25im*kout)*besselh.(1,kout*r).*((shape.ft[:,1].-ls_pos[1,1]).*shape.dft[:,2]-(shape.ft[:,2].-ls_pos[1,2]).*shape.dft[:,1])./r./ndft for i = 2:length(ls_amp) r = sqrt.((shape.ft[:,1] - ls_pos[i,1]).^2 + (shape.ft[:,2] - ls_pos[i,2]).^2) uls -= ls_amp[i]*0.25im*besselh.(0,kout*r) diff --git a/src/shapes.jl b/src/shapes.jl index cf4ab13..74b5487 100644 --- a/src/shapes.jl +++ b/src/shapes.jl @@ -6,7 +6,7 @@ Return a `ShapeParams` object containing the shape parametrized by """ function rounded_star(r, d, num, N) t = Float64[pi*j/N for j=0:(2*N-1)] - Rt = r + d*cos.(num*t) + Rt = r .+ d*cos.(num*t) dRt = (-d*num)*sin.(num*t) ft = [Rt.*cos.(t) Rt.*sin.(t)] dft = [dRt.*cos.(t) - Rt.*sin.(t) dRt.*sin.(t) + Rt.*cos.(t)] @@ -24,9 +24,9 @@ function squircle(r, N) t = Float64[pi*j/N for j=0:(2*N-1)] cost = cos.(t); cos4t = cos.(4*t); sint = sin.(t); sin4t = sin.(4*t); - rho = (sqrt(2)*r)./((3+cos4t).^0.25) + rho = (sqrt(2)*r)./((3 .+ cos4t).^0.25) ft = [rho.*cost rho.*sint] - drho = rho.*sin4t./(3+cos4t) + drho = rho.*sin4t./(3 .+ cos4t) dx = drho.*cost - rho.*sint dy = drho.*sint + rho.*cost dft = [dx dy] @@ -55,7 +55,7 @@ grid of points distanced `d`. function square_grid(a::Integer, d) offsetx = -0.5*(a-1) offsety = -0.5*(a-1) - centers = d*Float64[mod.((1:a^2)-1,a) + offsetx div.((1:a^2)-1,a) + offsety] + centers = d*[mod.(0:a^2-1, a) .+ offsetx div.(0:a^2-1, a) .+ offsety] end """ @@ -68,8 +68,8 @@ respectively. function rect_grid(a::Integer, b::Integer, dx, dy) offsetx = -0.5*(a-1) offsety = -0.5*(b-1) - xpoints = dx*((0:a-1) + offsetx) - ypoints = dy*((0:b-1) + offsety) + xpoints = dx*((0:a-1) .+ offsetx) + ypoints = dy*((0:b-1) .+ offsety) centers = [repeat(xpoints, outer=[b]) repeat(ypoints, inner=[a])] end @@ -84,7 +84,7 @@ If `minus1` is true, the last point in every odd row is omitted. function hex_grid(a::Integer, rows::Integer, d; minus1 = false) h = d*sqrt(0.75) #row height M = minus1 ? a*rows - div(rows,2) : a*rows - centers = Array{Float64}(M,2) + centers = Array{Float64}(undef, M, 2) ind = 1 for r = 0:rows-1 if mod(r,2) == 0 @@ -106,7 +106,7 @@ function hex_grid(a::Integer, rows::Integer, d; minus1 = false) end end end - offset = mean(centers, 1) + offset = mean(centers, dims=1) centers .-= offset end @@ -164,8 +164,8 @@ function randpoints(M, dmin, width, height, points; failures = 100) dist2 <= dmin2 && error("randpoints: given points have distance <= dmin") end - x_res = Array{Float64}(0) - y_res = Array{Float64}(0) + x_res = Array{Float64}(undef, 0) + y_res = Array{Float64}(undef, 0) fail = 0 while fail < failures && length(x_res) < M accepted = true @@ -260,14 +260,14 @@ is the radius of the rod centered at `(center[n,1],center[n,2])`. Otherwise quantizes the radii to uniformly spaced levels. """ function luneburg_grid(R_lens, N_cells, er; levels = 0, TM = true) - r_cell = Array{Float64}(N_cells^2) - centers = Array{Float64}(N_cells^2, 2) - flag_outside = Array{Bool}(N_cells^2) + r_cell = Array{Float64}(undef, N_cells^2) + centers = Array{Float64}(undef, N_cells^2, 2) + flag_outside = Array{Bool}(undef, N_cells^2) a = 2*R_lens/N_cells #cell dimension for ix = 1:N_cells, iy = 1:N_cells ind = (ix-1)*N_cells + iy - centers[ind,1] = a*(ix - N_cells/2.0 - 0.5); - centers[ind,2] = a*(iy - N_cells/2.0 - 0.5); + centers[ind,1] = a*(ix - N_cells/2.0 - 0.5) + centers[ind,2] = a*(iy - N_cells/2.0 - 0.5) rho2 = sum(abs2,centers[ind,:]) if rho2 > R_lens^2 #for now flag_outside[ind] = false @@ -281,7 +281,7 @@ function luneburg_grid(R_lens, N_cells, er; levels = 0, TM = true) if levels > 0 #uniform quantization - rs = collect(linspace(minimum(r_cell), maximum(r_cell), levels)) + rs = collect(range(minimum(r_cell), stop=maximum(r_cell), length=levels)) ddd = (r_cell - rs[1])/(rs[2] - rs[1]) ids = convert(Array{Int,1},round(ddd)) + 1 else diff --git a/src/utilities.jl b/src/utilities.jl index 4669871..ee4b959 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -63,12 +63,12 @@ well as a vector of indices `inds` of length `n` such that `v == u[inds]`. """ function uniqueind(v::Vector{T}) where T <: Number #returns inds,vals such that vals[inds[i]] == v[i] - inds = Array{Int}(0) - u = Array{T}(0) + inds = Array{Int}(undef, 0) + u = Array{T}(undef, 0) k = 0 for val in v - ind = findfirst(u, val) - if ind == 0 + ind = findfirst(isequal(val), u) + if ind == nothing k += 1 push!(u, val) push!(inds, k) @@ -86,8 +86,8 @@ Returns bounding box that contains all of the shapes in `sp`. """ function find_border(sp::ScatteringProblem) Rmax = maximum(s.R for s in sp.shapes) - x_max,y_max = maximum(sp.centers,1) + 2*Rmax - x_min,y_min = minimum(sp.centers,1) - 2*Rmax + x_max,y_max = maximum(sp.centers, dims=1) .+ 2*Rmax + x_min,y_min = minimum(sp.centers, dims=1) .- 2*Rmax border = [x_min; x_max; y_min; y_max] end @@ -99,8 +99,8 @@ Returns bounding box that contains all of the shapes in `sp` as well as specifie """ function find_border(sp::ScatteringProblem, points::Array{Float64,2}) Rmax = maximum(s.R for s in sp.shapes) - x_max,y_max = maximum([sp.centers;points],1) + 2*Rmax - x_min,y_min = minimum([sp.centers;points],1) - 2*Rmax + x_max,y_max = maximum([sp.centers;points], dims=1) .+ 2*Rmax + x_min,y_min = minimum([sp.centers;points], dims=1) .- 2*Rmax border = [x_min; x_max; y_min; y_max] end diff --git a/src/visualization.jl b/src/visualization.jl index 9420d1f..8c99f46 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -31,14 +31,14 @@ function plot_near_field(k0, kin, P, sp::ScatteringProblem, ui::Einc; # with pcolormesh, the field is calculated at the center of each rectangle. # thus we need two grids - the rectangle grid and the sampling grid. - x = linspace(x_min, x_max, x_points + 1) - y = linspace(y_min, y_max, y_points + 1) - xgrid = repmat(x', y_points + 1, 1) - ygrid = repmat(y, 1, x_points + 1) + x = range(x_min, stop=x_max, length=x_points + 1) + y = range(y_min, stop=y_max, length=y_points + 1) + xgrid = repeat(transpose(x), y_points + 1) + ygrid = repeat(y, 1, x_points + 1) dx = (x_max - x_min)/2/x_points dy = (y_max - y_min)/2/y_points - points = cat(2, vec(xgrid[1:y_points, 1:x_points]) + dx, - vec(ygrid[1:y_points, 1:x_points]) + dy) + points = cat(vec(xgrid[1:y_points, 1:x_points]) .+ dx, + vec(ygrid[1:y_points, 1:x_points]) .+ dy, dims=2) Ez = calc_near_field(k0, kin, P, sp, points, ui, method = method, opt = opt) @@ -51,7 +51,7 @@ function plot_near_field(k0, kin, P, sp::ScatteringProblem, ui::Einc; xlim([x_min/normalize;x_max/normalize]) ylim([y_min/normalize;y_max/normalize]) tight_layout() - ax[:set_aspect]("equal", adjustable = "box") + ax.set_aspect("equal", adjustable = "box") return (points,Ez),(xgrid,ygrid,zgrid) end @@ -69,13 +69,13 @@ function plot_far_field(k0, kin, P, sp::ScatteringProblem, pw::PlaneWave; Rmax = maximum(s.R for s in sp.shapes) - x_max,y_max = maximum(sp.centers,1) + 2*Rmax - x_min,y_min = minimum(sp.centers,1) - 2*Rmax + x_max,y_max = maximum(sp.centers, dims=1) + 2*Rmax + x_min,y_min = minimum(sp.centers, dims=1) - 2*Rmax Raggregate = 0.5*max(x_max - x_min, y_max - y_min) #radius of bounding circle x_center = 0.5*(x_max + x_min) y_center = 0.5*(y_max + y_min) Rfar = Raggregate*1e6 - theta_far = linspace(0, 2π, plot_points) + theta_far = range(0, stop=2π, length=plot_points) x_far = x_center + Rfar*cos.(theta_far) y_far = y_center + Rfar*sin.(theta_far) points = [x_far y_far] @@ -109,13 +109,13 @@ function draw_shapes(shapes, ids, centers, φs; ax = gca(), normalize = 1.0) if φs[ic] == 0.0 ft_rot = shapes[ids[ic]].ft .+ centers[ic,:]' else - Rot = cartesianrotation(φs[ic]) - ft_rot = shapes[ids[ic]].ft*Rot.' .+ centers[ic,:]' + Rot = cartesianrotation(φs[ic])' + ft_rot = shapes[ids[ic]].ft*Rot .+ centers[ic,:]' end - ax[:plot]([ft_rot[:,1];ft_rot[1,1]]/normalize, + ax.plot([ft_rot[:,1];ft_rot[1,1]]/normalize, [ft_rot[:,2];ft_rot[1,2]]/normalize, "k", linewidth = 2) else - ax[:add_patch](patch.Circle((centers[ic,1]/normalize, + ax.add_patch(patch.Circle((centers[ic,1]/normalize, centers[ic,2]/normalize), radius = shapes[ids[ic]].R/normalize, edgecolor="k", facecolor="none", linewidth = 2)) @@ -152,108 +152,109 @@ function calc_near_field(k0, kin, P, sp::ScatteringProblem, points, ui::Einc; result,sigma_mu = solve_particle_scattering_FMM(k0, kin, P, sp, ui, opt, verbose = verbose) if result[2].isconverged == false - warn("FMM process did not converge") + @warn("FMM process did not converge") return end beta = result[1] else - beta, sigma_mu = solve_particle_scattering(k0, kin, P, sp, ui) + beta, sigma_mu = solve_particle_scattering(k0, kin, P, sp, ui, + verbose = verbose) end - tic() #first, let's mark which points are in which shapes in tags: #0 denotes outside everything, +-i means inside shape i or between it and its "multipole disk" - tags = tagpoints(sp, points) - dt_tag = toq() + dt_tag = @elapsed begin + tags = tagpoints(sp, points) + end - tic() - rng_in = zeros(Bool,size(points,1)) - rng_out = zeros(Bool,size(points,1)) - Rot = Array{Float64}(2,2) - for ic = 1:size(sp) - rng_in[:] = (tags .== ic) - rng_out[:] = (tags .== -ic) - (any(rng_in) || any(rng_out)) || continue - if typeof(shapes[ids[ic]]) == ShapeParams - if φs[ic] == 0.0 - ft_rot = shapes[ids[ic]].ft .+ centers[ic,:]' - dft_rot = shapes[ids[ic]].dft - else - Rot[:] = cartesianrotation(φs[ic]) - ft_rot = shapes[ids[ic]].ft*Rot.' .+ centers[ic,:]' - dft_rot = shapes[ids[ic]].dft*Rot.' - end - end - #field inside shape - if any(rng_in) + dt_in = @elapsed begin + rng_in = zeros(Bool,size(points,1)) + rng_out = zeros(Bool,size(points,1)) + Rot = Array{Float64}(undef, 2,2) + for ic = 1:size(sp) + rng_in[:] = (tags .== ic) + rng_out[:] = (tags .== -ic) + (any(rng_in) || any(rng_out)) || continue if typeof(shapes[ids[ic]]) == ShapeParams - u[rng_in] += scatteredfield(sigma_mu[ic], kin, shapes[ids[ic]].t, - ft_rot, dft_rot, points[rng_in,:]) - else - u[rng_in] += innerFieldCircle(kin, sigma_mu[ic], centers[ic,:], - points[rng_in,:]) + if φs[ic] == 0.0 + ft_rot = shapes[ids[ic]].ft .+ centers[ic,:]' + dft_rot = shapes[ids[ic]].dft + else + Rot[:] = cartesianrotation(φs[ic])' + ft_rot = shapes[ids[ic]].ft*Rot .+ centers[ic,:]' + dft_rot = shapes[ids[ic]].dft*Rot + end end - end - #field between shape and multipole disk (impossible for circle) - if any(rng_out) - u[rng_out] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t, - ft_rot, dft_rot, points[rng_out,:]) - u[rng_out] += uinc(k0, points[rng_out,:], ui) - for ic2 = 1:size(sp) - ic == ic2 && continue - if method == "multipole" - scattered_field_multipole!(u, k0, beta, P, centers, ic2, points, - find(rng_out)) - elseif method == "recurrence" - scattered_field_multipole_recurrence!(u, k0, beta, P, centers, ic2, points, - find(rng_out)) + #field inside shape + if any(rng_in) + if typeof(shapes[ids[ic]]) == ShapeParams + u[rng_in] += scatteredfield(sigma_mu[ic], kin, shapes[ids[ic]].t, + ft_rot, dft_rot, points[rng_in,:]) else - if φs[ic2] == 0.0 - ft_rot2 = shapes[ids[ic2]].ft .+ centers[ic2,:]' - dft_rot2 = shapes[ids[ic2]].dft + u[rng_in] += innerFieldCircle(kin, sigma_mu[ic], centers[ic,:], + points[rng_in,:]) + end + end + #field between shape and multipole disk (impossible for circle) + if any(rng_out) + u[rng_out] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t, + ft_rot, dft_rot, points[rng_out,:]) + u[rng_out] += uinc(k0, points[rng_out,:], ui) + for ic2 = 1:size(sp) + ic == ic2 && continue + if method == "multipole" + scattered_field_multipole!(u, k0, beta, P, centers, ic2, points, + findall(rng_out)) + elseif method == "recurrence" + scattered_field_multipole_recurrence!(u, k0, beta, P, centers, ic2, points, + findall(rng_out)) else - Rot[:] = cartesianrotation(φs[ic2]) - ft_rot2 = shapes[ids[ic2]].ft*Rot.' .+ centers[ic2,:]' - dft_rot2 = shapes[ids[ic2]].dft*Rot.' + if φs[ic2] == 0.0 + ft_rot2 = shapes[ids[ic2]].ft .+ centers[ic2,:]' + dft_rot2 = shapes[ids[ic2]].dft + else + Rot[:] = cartesianrotation(φs[ic2])' + ft_rot2 = shapes[ids[ic2]].ft*Rot .+ centers[ic2,:]' + dft_rot2 = shapes[ids[ic2]].dft*Rot + end + u[rng_out] += scatteredfield(sigma_mu[ic2], k0, shapes[ids[ic2]].t, + ft_rot2, dft_rot2, points[rng_out,:]) end - u[rng_out] += scatteredfield(sigma_mu[ic2], k0, shapes[ids[ic2]].t, - ft_rot2, dft_rot2, points[rng_out,:]) end end end end - dt_in = toq() - tic() #now compute field outside all shapes - rng = (tags .== 0) - #incident field - u[rng] = uinc(k0, points[rng,:], ui) - if method == "multipole" - scattered_field_multipole!(u, k0, beta, P, centers, 1:size(sp), points, find(rng)) - elseif method == "recurrence" - scattered_field_multipole_recurrence!(u, k0, beta, P, centers, 1:size(sp), points, find(rng)) - else - for ic = 1:size(centers,1) - if typeof(shapes[ids[ic]]) == ShapeParams - if φs[ic] == 0.0 - ft_rot = shapes[ids[ic]].ft .+ centers[ic,:]' - u[rng] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t, - ft_rot, shapes[ids[ic]].dft, points[rng,:]) + dt_out = @elapsed begin + rng = (tags .== 0) + #incident field + u[rng] = uinc(k0, points[rng,:], ui) + if method == "multipole" + scattered_field_multipole!(u, k0, beta, P, centers, 1:size(sp), points, findall(rng)) + elseif method == "recurrence" + scattered_field_multipole_recurrence!(u, k0, beta, P, centers, 1:size(sp), points, findall(rng)) + else + for ic = 1:size(centers,1) + if typeof(shapes[ids[ic]]) == ShapeParams + if φs[ic] == 0.0 + ft_rot = shapes[ids[ic]].ft .+ centers[ic,:]' + u[rng] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t, + ft_rot, shapes[ids[ic]].dft, points[rng,:]) + else + Rot[:] = copy(transpose(cartesianrotation(φs[ic]))) + ft_rot = shapes[ids[ic]].ft*Rot .+ centers[ic,:]' + dft_rot = shapes[ids[ic]].dft*Rot + u[rng] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t, + ft_rot, dft_rot, points[rng,:]) + end else - Rot[:] = cartesianrotation(φs[ic]) - ft_rot = shapes[ids[ic]].ft*Rot.' .+ centers[ic,:]' - dft_rot = shapes[ids[ic]].dft*Rot.' - u[rng] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t, - ft_rot, dft_rot, points[rng,:]) + scattered_field_multipole!(u, k0, beta, P, centers, ic, points, + findall(rng)) end - else - scattered_field_multipole!(u, k0, beta, P, centers, ic, points, - find(rng)) end end end - dt_out = toq() if verbose println("Time spent calculating field:") println("Location tagging: $dt_tag") @@ -270,7 +271,7 @@ function calc_far_field(k0, kin, P, points, sp::ScatteringProblem, pw::PlaneWave if opt.FMM result,sigma_mu = solve_particle_scattering_FMM(k0, kin, P, sp, pw, opt) if result[2].isconverged == false - warn("FMM process did not converge") + @warn("FMM process did not converge") return end beta = result[1] @@ -293,9 +294,9 @@ function calc_far_field(k0, kin, P, points, sp::ScatteringProblem, pw::PlaneWave Ez[:] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t, ft_rot, shapes[ids[ic]].dft, points) else - Rot = cartesianrotation(φs[ic]) - ft_rot = shapes[ids[ic]].ft*Rot.' .+ centers[ic,:]' - dft_rot = shapes[ids[ic]].dft*Rot.' + Rot = copy(transpose(cartesianrotation(φs[ic]))) + ft_rot = shapes[ids[ic]].ft*Rot .+ centers[ic,:]' + dft_rot = shapes[ids[ic]].dft*Rot Ez[:] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t, ft_rot, dft_rot, points) end @@ -313,7 +314,7 @@ function tagpoints_old(sp, points) shapes = sp.shapes; ids = sp.ids; centers = sp.centers; φs = sp.φs tags = zeros(Integer, size(points,1)) - X = Array{Float64}(2) + X = Array{Float64}(undef, 2) for ix = 1:size(points,1) for ic = 1:size(sp) X .= points[ix,:] - centers[ic,:] @@ -339,7 +340,7 @@ function tagpoints(sp, points) shapes = sp.shapes; ids = sp.ids; centers = sp.centers; φs = sp.φs tags = zeros(Integer, size(points,1)) - X = Array{Float64}(2) #tmp arrays + X = Array{Float64}(undef, 2) #tmp arrays for ix = 1:size(points,1) #need two loops due to "break" for ic = 1:size(sp) X[1] = points[ix,1] - centers[ic,1] diff --git a/test/fmm_test.jl b/test/fmm_test.jl index 4771c9d..b4cd2f4 100644 --- a/test/fmm_test.jl +++ b/test/fmm_test.jl @@ -14,19 +14,19 @@ φs = 2π*rand(M^2) #random rotation angles dist = 2*maximum(shapes[i].R for i=1:2) try - centers = randpoints(M^2, dist, 5λ0, 5λ0, [0.0 0.0; 0.01+dist 0.01], - failures = 1_000) - sp = ScatteringProblem(shapes, ids, centers, φs) - @test verify_min_distance(sp, [0.0 0.0; 0.01+dist 0.01]) + centers = randpoints(M^2, dist, 5λ0, 5λ0, [0.0 0.0; 0.01+dist 0.01], + failures = 1_000) + sp = ScatteringProblem(shapes, ids, centers, φs) + @test verify_min_distance(sp, [0.0 0.0; 0.01+dist 0.01]) catch - warn("Could not find random points (1)") + @warn("Could not find random points (1)") end try - centers = randpoints(M^2, dist, 5λ0, 5λ0, failures = 1_000) - sp = ScatteringProblem(shapes, ids, centers, φs) - @test verify_min_distance(sp) + centers = randpoints(M^2, dist, 5λ0, 5λ0, failures = 1_000) + sp = ScatteringProblem(shapes, ids, centers, φs) + @test verify_min_distance(sp) catch - warn("Could not find random points (2)") + @warn("Could not find random points (2)") end centers = rect_grid(M, M, 0.8λ0, λ0) @@ -39,16 +39,31 @@ mFMM2 = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize; tri = true) @test norm(mFMM1.Znear - mFMM2.Znear, Inf) == 0 - points = [linspace(-λ0*M, λ0*M, 200) zeros(200)] + points = [range(-λ0*M, stop=λ0*M, length=200) zeros(200)] u1 = calc_near_field(k0, kin, P, sp, points, PlaneWave(θ_i); opt = fmm_options, verbose = false) u2 = calc_near_field(k0, kin, P, sp, points, PlaneWave(θ_i); opt = fmm_options, verbose = true, method = "density") @test norm(u1 - u2)/norm(u1) < 1e-6 - u3 = calc_near_field(k0, kin, P, sp, points, LineSource(-0.8λ0,-λ0); opt = fmm_options, - verbose = false) - u4 = calc_near_field(k0, kin, P, sp, points, LineSource(-0.8λ0,-λ0); opt = fmm_options, - verbose = true, method = "recurrence") + u3 = calc_near_field(k0, kin, P, sp, points, LineSource(-0.8λ0,-λ0); + opt = fmm_options, verbose = false) + u4 = calc_near_field(k0, kin, P, sp, points, LineSource(-0.8λ0,-λ0); + opt = fmm_options, verbose = true, method = "recurrence") @test norm(u3 - u4)/norm(u3) < 1e-6 + + #test transpose fmm + mFMM, sMatrices, sLU, buf = ParticleScattering.prepare_fmm_reusal_φs(k0, + kin, P, shapes, centers, ids, fmm_options) + v1 = 5*rand(ComplexF64, M^2*(2P+1)) + v2 = rand(ComplexF64, M^2*(2P+1)) + o1 = Array{ComplexF64}(undef, M^2*(2P+1)) + o2 = Array{ComplexF64}(undef, M^2*(2P+1)) + ParticleScattering.FMM_mainMVP_pre!(o1, v1, sMatrices, φs, ids, P, mFMM, + buf.pre_agg, buf.trans) + res1 = transpose(v2)*o1 + ParticleScattering.FMM_mainMVP_transpose!(o2, v2, sMatrices, φs, ids, P, + mFMM, buf.pre_agg, buf.trans) + res2 = transpose(v1)*o2 + @test res1 ≈ res2 end diff --git a/test/multipole_test.jl b/test/multipole_test.jl index 27779e9..2cf4b1f 100644 --- a/test/multipole_test.jl +++ b/test/multipole_test.jl @@ -7,7 +7,7 @@ N = 1000 P = 5 - shapes = [squircle(2π/λ0, N)] + shapes = [squircle(0.2λ0, N)] centers = [-1.5λ0 0; 1.5λ0 0] Ns = size(centers,1) φs = [0.0; 0.23] @@ -23,11 +23,11 @@ #now rotate everything by some angle and compare θ_r = 1.23#2π*rand() #notation is transposed due to structure of centers - shapes2 = [squircle(2π/λ0, N); CircleParams(15)] #for extra code coverage - centers2 = centers*([cos(θ_r) -sin(θ_r); sin(θ_r) cos(θ_r)].') - φs2 = θ_r + φs + shapes2 = [squircle(0.2λ0, N); CircleParams(15)] #no effect,for extra code coverage + centers2 = centers*[cos(θ_r) sin(θ_r); -sin(θ_r) cos(θ_r)] + φs2 = θ_r .+ φs sp2 = ScatteringProblem(shapes2, ids, centers2, φs2) - points2 = points*([cos(θ_r) -sin(θ_r); sin(θ_r) cos(θ_r)].') + points2 = points*[cos(θ_r) sin(θ_r); -sin(θ_r) cos(θ_r)] β2, σ2 = solve_particle_scattering(k0, kin, P, sp2, PlaneWave(θ_i + θ_r); get_inner = true, verbose = true) diff --git a/test/optimize_angle_pwr_test.jl b/test/optimize_angle_pwr_test.jl new file mode 100644 index 0000000..1020e54 --- /dev/null +++ b/test/optimize_angle_pwr_test.jl @@ -0,0 +1,61 @@ +@testset "optimize angle power" begin + λ0 = 1 #doesn't matter since everything is normalized to λ0 + k0 = 2π/λ0 + kin = 3k0 + ui = PlaneWave(0.0) + M = 10 + + sfun(N) = rounded_star(0.35λ0, 0.1λ0, 4, N) + # N, errN = minimumN(k0, kin, sfun; tol=1e-6, N_start = 200, N_min = 100, + # N_max = 400) + shapes = [sfun(202)] + P = 12#P, errP = minimumP(k0, kin, shapes[1]; P_max = 30, tol = 1e-7) + + centers = rect_grid(2, div(M,2), λ0, λ0) #2xM/2 grid with distance λ0 + ids = ones(Int, M) + φs0 = zeros(M) + sp = ScatteringProblem(shapes, ids, centers, φs0) + + points = [2λ0*ones(15) range(minimum(centers[:,2]), stop=maximum(centers[:,2]), length=15)] + nhat = [ones(15) zeros(15)] + len = norm(points[1,:] - points[end,:]) + @assert verify_min_distance(sp, points) + + fmm_options = FMMoptions(true, acc = 6, dx = 2λ0) + groups,boxSize = divideSpace(sp.centers, fmm_options) + P2,Q = FMMtruncation(fmm_options.acc, boxSize, k0) + optim_options = Optim.Options(f_tol = 1e-5, iterations = 50, + store_trace = true, show_trace = false) + + # Allocate buffers + opb = [OptimProblemBuffer(k0, kin, sp.centers, ui, P)] + power_buffer = [optMatrixPwr(points, sp.centers, M, P, k0, ui, 1, nhat, len)] + mFMM = [FMMbuildMatrices(k0, P, P2, Q, groups, sp.centers, boxSize, tri = true)] + buf = FMMbuffer(M, P, Q, length(groups)) + scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids) + scatteringMatrices = [scatteringMatrices] + scatteringLU = [[lu(scatteringMatrices[1][iid]) for iid = 1:length(shapes)]] + + function fobj_test(sv::Array{PowerBuffer}) + res = -sv[1].pow + end + + function gobj_test!(sv::Array{PowerBuffer}, opb::Array{OptimProblemBuffer}) + #calculate -(∂f/∂β)ᵀ for adjoint method + opb[1].rhs_grad[:] = sv[1].∂pow + end + + #initial, last must be different before first iteration + last_φs = similar(φs0); last_φs[1] = NaN; @assert(last_φs != φs0) + df = Optim.OnceDifferentiable( + φs -> optimize_pwr_φ_f(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, M, buf, fmm_options, mFMM, fobj_test), + (grad_stor, φs) -> optimize_pwr_φ_g!(grad_stor, φs, last_φs, opb, power_buffer, ids, scatteringMatrices, scatteringLU, P, M, buf, fmm_options, mFMM, gobj_test!), + φs0) + + res = Optim.optimize(df, φs0, Optim.LBFGS(), optim_options) + + power_after = calc_power(k0, kin, P, ScatteringProblem(shapes, ids, centers, + res.minimizer), points, nhat, ui)*len + @test res.f_converged + @test power_after ≈ 0.00094388385 +end diff --git a/test/optimize_angle_test.jl b/test/optimize_angle_test.jl index 49340d3..0895aae 100644 --- a/test/optimize_angle_test.jl +++ b/test/optimize_angle_test.jl @@ -1,6 +1,4 @@ #based on the optimizing angle tutorial -import Optim - @testset "optimize angle" begin λ0 = 1 #doesn't matter since everything is normalized to λ0 k0 = 2π/λ0 @@ -20,7 +18,7 @@ import Optim sp = ScatteringProblem(shapes, ids, centers, φs0) points = [-0.05 0.0; 0.05 0.0; 0.0 0.05; 0.0 -0.05] - assert(verify_min_distance(sp, points)) + @assert verify_min_distance(sp, points) fmm_options = FMMoptions(true, acc = 6, dx = 2λ0) diff --git a/test/optimize_radius_pwr_test.jl b/test/optimize_radius_pwr_test.jl new file mode 100644 index 0000000..c640cd2 --- /dev/null +++ b/test/optimize_radius_pwr_test.jl @@ -0,0 +1,83 @@ +#two-sided diode-type test. +@testset "optimize radius power" begin + λ0 = 1 #doesn't matter since everything is normalized to λ0 + k0 = 2π/λ0 + kin = 4k0 + ui = [PlaneWave(); PlaneWave(π)] + a = λ0/2 + P = 6 + + centers = rect_grid(4, 8, a, a) #symmetric 4x8 grid with distance a + M = size(centers, 1) + ids = ones(Int, M) + φs = zeros(M) + centers_abs = centers[:,1] + 1im*abs.(centers[:,2]) + ids, centers_abs = ParticleScattering.uniqueind(centers_abs) + Nrs = maximum(ids) #if some are the same, Nrs is the number of different r + initial_rs = 0.25*a*ones(Nrs) + r_max = 0.45*a*ones(Nrs) + r_min = 1e-4*a*ones(Nrs) + shapes = CircleParams.(initial_rs) + sp = ScatteringProblem(shapes, ids, centers, φs) + Npoints = 15 + points1 = [2λ0*ones(Npoints) range(-λ0, stop=λ0, length=Npoints)] + points2 = [-2λ0*ones(Npoints) range(-λ0, stop=λ0, length=Npoints)] + nhat1 = [ones(Npoints) zeros(Npoints)] + nhat2 = [-ones(Npoints) zeros(Npoints)] + len = norm(points1[1,:] - points1[end,:]) #same for 2 + @assert verify_min_distance(sp, [points1;points2]) + + fmm_options = FMMoptions(true, acc = 6, dx = λ0) + groups,boxSize = divideSpace(sp.centers, fmm_options) + P2,Q = FMMtruncation(fmm_options.acc, boxSize, k0) + optim_options = Optim.Options(f_tol = 1e-5, iterations = 10, outer_iterations = 5, + store_trace = false, show_trace = false) + # Allocate buffers + opb = [OptimProblemBuffer(k0, kin, sp.centers, ui[1], P), + OptimProblemBuffer(k0, kin, sp.centers, ui[2], P)] + power_buffer = [optMatrixPwr(points1, sp.centers, M, P, k0, ui[1], 1, nhat1, len), + optMatrixPwr(points2, sp.centers, M, P, k0, ui[2], 2, nhat2, len)] + mFMM = [FMMbuildMatrices(k0, P, P2, Q, groups, sp.centers, boxSize, tri = true) for i=1:2] + buf = FMMbuffer(M, P, Q, length(groups)) + scatteringMatrices = [[sparse(one(Complex{Float64})I, 2P+1, 2P+1) for i = 1:Nrs] for i = 1:2] + dS_S = [[sparse(one(Complex{Float64})I, 2P+1, 2P+1) for i = 1:Nrs] for i = 1:2] + + function fobj_testr(sv::Array{PowerBuffer}) + if sv[1].pow < 0 || sv[2].pow < 0 + return abs(sv[2].pow/sv[1].pow)*1e9 #instead of Inf + end + barr = -(log(sv[1].pow) + log(sv[2].pow)) #so they're > 0 + res = sv[2].pow/sv[1].pow + res*(1 + barr) + end + + function gobj_testr!(sv::Array{PowerBuffer}, opb::Array{OptimProblemBuffer}) + #calculate -(∂f/∂β)ᵀ for adjoint method + if sv[1].pow < 0 || sv[2].pow < 0 + opb[1].rhs_grad[:] .= 0 + opb[2].rhs_grad[:] .= 0 + return + end + barr = -(log(sv[1].pow) + log(sv[2].pow)) + res = sv[2].pow/sv[1].pow + opb[1].rhs_grad[:] = sv[1].∂pow*((1 + barr)*sv[2].pow/sv[1].pow^2) + opb[2].rhs_grad[:] = sv[2].∂pow*(-(1 + barr)/sv[1].pow) + opb[1].rhs_grad .+= sv[1].∂pow*(res/sv[1].pow) + opb[2].rhs_grad .+= sv[2].∂pow*(res/sv[2].pow) + end + + #initial, last must be different before first iteration + last_rs = similar(initial_rs); @assert last_rs != initial_rs #initial_rs, last_rs must be different before first iteration + df = Optim.OnceDifferentiable( + rs -> optimize_pwr_rs_f(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, M, buf, fmm_options, φs, mFMM, fobj_testr), + (grad_stor, rs) -> optimize_pwr_rs_g!(grad_stor, rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, M, buf, fmm_options, φs, mFMM, gobj_testr!), + initial_rs) + + res = Optim.optimize(df, r_min, r_max, initial_rs, Optim.Fminbox(Optim.LBFGS()), optim_options) + sp_after = ScatteringProblem(CircleParams.(res.minimizer), ids, centers, φs) + power_after_right = calc_power(k0, kin, P, sp_after, points1, nhat1, ui[1])*len + power_after_left = calc_power(k0, kin, P, sp_after, points2, nhat2, ui[2])*len + @test power_after_right > 0 + @test power_after_left > 0 + @test power_after_right/power_after_left > 1 +end diff --git a/test/optimize_radius_test.jl b/test/optimize_radius_test.jl index 830ea80..35ceda7 100644 --- a/test/optimize_radius_test.jl +++ b/test/optimize_radius_test.jl @@ -1,6 +1,4 @@ #based on the optimizing radius tutorial -import Optim - @testset "optimize radius" begin er = 4.5 k0 = 2π @@ -26,8 +24,8 @@ import Optim r_max = (0.4*a)*ones(J) r_min = (1e-3*a)*ones(J) rs0 = (0.25*a)*ones(J) - assert(verify_min_distance([CircleParams(r_max[i]) for i = 1:J], - centers, ids, points)) + @assert verify_min_distance([CircleParams(r_max[i]) for i = 1:J], + centers, ids, points) res = optimize_radius(rs0, r_min, r_max, points, ids, P, ui, k0, kin, centers, fmm_options, optim_options, minimize = true) diff --git a/test/runtests.jl b/test/runtests.jl index 70803cd..3860691 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,10 +1,13 @@ using ParticleScattering -using Base.Test +using Test, LinearAlgebra +import Optim, LineSearches, SparseArrays.sparse include("scatteredfield_test.jl") include("multipole_test.jl") include("fmm_test.jl") include("optimize_radius_test.jl") include("optimize_angle_test.jl") +include("optimize_angle_pwr_test.jl") +include("optimize_radius_pwr_test.jl") include("minimum_N_P_test.jl") include("utilities_test.jl") diff --git a/test/utilities_test.jl b/test/utilities_test.jl index 4794ca5..598d59b 100644 --- a/test/utilities_test.jl +++ b/test/utilities_test.jl @@ -2,7 +2,7 @@ s = rounded_star(1, 0.1, 7, 100) @test ParticleScattering.pInPolygon([0.0 0.0], s.ft) @test !ParticleScattering.pInPolygon([0.84 0.35], s.ft) - @test ParticleScattering.pInPolygon([0.84 0.35], s.ft+0.1) + @test ParticleScattering.pInPolygon([0.84 0.35], s.ft .+ 0.1) sp = ScatteringProblem([s],[1],[0.1 0.1], [15π]) border = find_border(sp) @@ -10,4 +10,9 @@ @test border == border2 border3 = find_border(sp,[0.1 0.9999]) @test border != border3 + + vs = [1.1im, -1.2, 1 + 324.0im, 1.1im, 1 - 324im, 1.0 + 324im] + is, vs2 = uniqueind(vs) + @test all(vs2[is] .== vs) + @test length(vs2) == 4 end