diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..3b3a0f1
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*.dat
+.ipynb_checkpoints/
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..94a9ed0
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,674 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc.
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+
+ Copyright (C)
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see .
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+ Copyright (C)
+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+.
+
+ The GNU General Public License does not permit incorporating your program
+into proprietary programs. If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License. But first, please read
+.
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..e3c14a0
--- /dev/null
+++ b/README.md
@@ -0,0 +1,36 @@
+# HelmholtzMOR
+Module for the solution and model order reduction of the Helmholtz problem. Coded in Python 3.6 and based on Fenics.
+
+## Installing
+After cloning the repository, you will need to have Fenics installed on your system. To this aim, you can install Anaconda3/Miniconda3 and run the command (from the main directory ./)
+```
+conda env create --file conda-fenics-nonotebook.yml
+```
+or
+```
+conda env create --file conda-fenics.yml
+```
+
+This will create an environment where Fenics can be used. To activate the environment, it is sufficent to run (from any directory)
+```
+source activate fenicsenv-nonotebook
+```
+or
+```
+source activate fenicsenv
+```
+To deactivate it, run
+```
+source deactivate
+```
+
+## Running the samples
+Many examples can be found in the ./examples/ folder. Both Python and IPhyton scripts are given.
+To run the Python examples you can use
+```
+python examplename.py
+```
+whereas IPython examples can be run (e.g.) through Jupyter notebook. In both cases, an environment/kernel with Fenics installed must be active.
+
+## License
+This project is licensed under the GNU GENERAL PUBLIC LICENSE license - see the [LICENSE.md](LICENSE.md) file for details.
diff --git a/examples/Python/HelmholtzLagrangeApproximantsSweep.py b/examples/Python/HelmholtzLagrangeApproximantsSweep.py
new file mode 100644
index 0000000..2a9d347
--- /dev/null
+++ b/examples/Python/HelmholtzLagrangeApproximantsSweep.py
@@ -0,0 +1,100 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Example homogeneous Dirichlet forcing wave SWEEP
+
+from __future__ import print_function
+import fenics as fen
+import numpy as np
+import sympy as sp
+from context import FenicsHelmholtzEngine as HFEngine
+from context import FenicsHSEngine as HSEngine
+from context import ROMApproximantLagrangePade as Pade
+from context import ROMApproximantLagrangeRB as RB
+from context import ROMApproximantSweeper as Sweeper
+
+PI = np.pi
+
+nu = 12**.5
+theta = PI / 3
+npoints = 31
+ktars = np.linspace(0, 21, npoints)
+
+x, y = sp.symbols('x[0] x[1]', real=True)
+wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
+phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
+uex = wex * sp.exp(-1.j * phiex)
+fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
+
+nx = ny = 10
+mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+solver = HFEngine.FenicsHelmholtzEngine(mesh = mesh, wavenumber = nu, forcingTerm = forcingTerm, FEDegree = 3,
+ DirichletBoundary = 'all', DirichletDatum = 0)
+plotter = HSEngine.FenicsHSEngine(solver.V)
+
+shift = 3
+nsets = 5
+stride = 3
+Smax = stride * (nsets - 1) + shift + 2
+paramsPade = {'S':Smax, 'polyBasis':'CHEBYSHEV', 'POD':True}
+paramsRB = {'S':Smax, 'nodesType':'CHEBYSHEV', 'POD':True}
+paramsSetsPade = [None] * nsets
+paramsSetsRB = [None] * nsets
+for i in range(nsets):
+ paramsSetsPade[i] = {'N': stride * i + shift + 1, 'M': stride * i + shift,
+ 'S': stride * i + shift + 2}
+ paramsSetsRB[i] = {'R': stride * i + shift + 1, 'S': stride * i + shift + 2}
+
+appPade = Pade.ROMApproximantLagrangePade(solver, plotter, ks = [4 + .5j, 14 + .5j],
+ w = nu, approxParameters = paramsPade)
+appRB = RB.ROMApproximantLagrangeRB(solver, plotter, ks = [4 + .5j, 14 + .5j],
+ w = nu, approxParameters = paramsRB)
+
+sweeper = Sweeper.ROMApproximantSweeper(ktars = ktars, mostExpensive = 'Approx')
+sweeper.ROMEngine = appPade
+sweeper.params = paramsSetsPade
+filenamePade = sweeper.sweep('../Data/HelmholtzBubbleLagrangePade.dat', outputs = 'ALL')
+
+sweeper.ROMEngine = appRB
+sweeper.params = paramsSetsRB
+filenameRB = sweeper.sweep('../Data/HelmholtzBubbleLagrangeRB.dat', outputs = 'ALL')
+
+####################
+
+from matplotlib import pyplot as plt
+plt.jet()
+
+for i in range(nsets):
+ nSamples = stride*i+shift+2
+ PadeOutput = sweeper.read(filenamePade, {'S':nSamples},
+ ['kRe', 'HFNorm', 'AppNorm', 'AppError'])
+ RBOutput = sweeper.read(filenameRB, {'S':nSamples},
+ ['kRe', 'AppNorm', 'AppError'])
+
+ ktarsF = PadeOutput['kRe']
+ solNormF = PadeOutput['HFNorm']
+ PadektarsF = PadeOutput['kRe']
+ PadeNormF = PadeOutput['AppNorm']
+ PadeErrorF = PadeOutput['AppError']
+ RBktarsF = RBOutput['kRe']
+ RBNormF = RBOutput['AppNorm']
+ RBErrorF = RBOutput['AppError']
+
+ plt.figure()
+ plt.semilogy(ktarsF, solNormF, 'k-', label='Sol norm')
+ plt.semilogy(PadektarsF, PadeNormF, 'b.--', label='Pade'' norm, S = {}'.format(nSamples))
+ plt.semilogy(RBktarsF, RBNormF, 'g.--', label='RB norm, S = {}'.format(nSamples))
+ plt.legend()
+ plt.grid()
+ plt.figure()
+ plt.semilogy(PadektarsF, PadeErrorF, 'b', label='Pade'' error, S = {}'.format(nSamples))
+ plt.semilogy(RBktarsF, RBErrorF, 'g', label='RB error, S = {}'.format(nSamples))
+ plt.legend()
+ plt.grid()
+ plt.figure()
+ plt.semilogy(ktarsF, PadeErrorF / solNormF, 'b', label='Pade'' relative error, S = {}'.format(nSamples))
+ plt.semilogy(RBktarsF, RBErrorF / solNormF, 'g', label='RB relative error, S = {}'.format(nSamples))
+ plt.legend()
+ plt.grid()
diff --git a/examples/Python/HelmholtzPadeLagrangeApproximant.py b/examples/Python/HelmholtzPadeLagrangeApproximant.py
new file mode 100644
index 0000000..8664d06
--- /dev/null
+++ b/examples/Python/HelmholtzPadeLagrangeApproximant.py
@@ -0,0 +1,181 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from __future__ import print_function
+import fenics as fen
+import numpy as np
+import sympy as sp
+from context import FenicsHelmholtzEngine as HFEngine
+from context import FenicsHSEngine as HSEngine
+from context import ROMApproximantLagrangePade as Pade
+
+PI = np.pi
+
+testNo = 4
+
+if testNo == 1:
+
+ params = {'N':4, 'M':3, 'S':5, 'polyBasis':'CHEBYSHEV', 'POD':True}
+
+ nu = 12**.5
+ theta = PI / 3
+ ztar = 11
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
+ phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
+ uex = wex * sp.exp(-1.j * phiex)
+ fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
+
+ nx = ny = 40
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+ forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+ solver = HFEngine.FenicsHelmholtzEngine(mesh = mesh, wavenumber = nu, forcingTerm = forcingTerm, FEDegree = 3,
+ DirichletBoundary = 'all', DirichletDatum = 0)
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = Pade.ROMApproximantLagrangePade(solver, plotter, ks = [10 + .5j, 14 + .5j],
+ w = np.real(nu), approxParameters = params)
+
+ approx.plotApp(ztar, name = 'u_Pade''')
+ approx.plotHF(ztar, name = 'u_HF')
+ approx.plotErr(ztar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
+ print(appErr, solNorm, appErr/solNorm)
+
+ print(approx.getPoles(True))
+
+############
+elif testNo == 2:
+
+ params = {'N':9, 'M':8, 'S':10, 'polyBasis':'CHEBYSHEV', 'POD':True}
+ ztar = 3.9 ** 2.
+
+ n1 = 2**.5
+ n2 = 3**.5
+ kappa = 4
+ theta = PI * 75 / 180
+ d1, d2 = np.cos(theta), np.sin(theta)
+ K1 = kappa * n1 * d1
+ if kappa * n2 >= K1:
+ K2 = ((kappa*n2)**2 - K1**2)**.5
+ else:
+ K2 = 1.j * (K1**2 - (kappa*n2)**2)**.5
+ R = (kappa * n1 * d2 - K2) / (kappa * n1 * d2 + K2)
+ T = R + 1
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ uex1 = T*sp.exp(1.j*(K1*x+K2*y))
+ uex2 = sp.exp(1.j*kappa*n1*(d1*x+d2*y)) + R*sp.exp(1.j*kappa*n1*(d1*x-d2*y))
+
+ # Exact solution
+ uexRe = fen.Expression('x[1]>=0 ? {0} : {1}'.format(\
+ sp.printing.ccode(sp.re(uex1)), sp.printing.ccode(sp.re(uex2))), degree=4)
+ uexIm = fen.Expression('x[1]>=0 ? {0} : {1}'.format(\
+ sp.printing.ccode(sp.im(uex1)), sp.printing.ccode(sp.im(uex2))), degree=4)
+
+ # wavenumber term
+ nRe = fen.Expression('x[1]<0 ? n1r : n2r', n1r = n1.real, n2r = n2.real, degree=4)
+ nIm = fen.Expression('x[1]<0 ? n1i : n2i', n1i = n1.imag, n2i = n2.imag, degree=4)
+
+ # Create mesh and define function space
+ nx = ny = 50
+ mesh = fen.RectangleMesh(fen.Point(-PI/2,-PI/2), fen.Point(PI/2,PI/2), nx, ny)
+
+ solver = HFEngine.FenicsHelmholtzEngine(mesh = mesh, wavenumber = kappa, refractionIndex = (nRe, nIm),
+ forcingTerm = 0, FEDegree = 3, DirichletBoundary = 'all',
+ DirichletDatum = (uexRe, uexIm))
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = Pade.ROMApproximantLagrangePade(solver, plotter, ks = np.power([3.85 + .15j, 4.15 + .15j], 2.),
+ w = kappa, approxParameters = params, plotSnapshots = True)
+
+ approx.plotApp(ztar, name = 'u_Pade''')
+ approx.plotHF(ztar, name = 'u_HF')
+ approx.plotErr(ztar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
+ print(appErr, solNorm, appErr/solNorm)
+
+ print(approx.getPoles(True))
+
+############
+elif testNo == 3:
+
+ def Dboundary(x, on_boundary):
+ return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
+
+ A = 10
+ kappa = 4
+ theta = PI * 90 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ u0ex = - A * sp.exp(1.j * phiex)
+
+ nx = ny = 40
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+
+ DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
+
+ params = {'N':40, 'M':39, 'S':45, 'polyBasis':'CHEBYSHEV', 'POD':True}
+
+ solver = HFEngine.FenicsHelmholtzScatteringEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0,
+ FEDegree = 3, DirichletBoundary = Dboundary,
+ RobinBoundary = 'rest', DirichletDatum = DirichletTerm)
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = Pade.ROMApproximantLagrangePade(solver, plotter, ks = [0, 8],
+ approxParameters = params,
+ plotSnapshots = False)
+
+ print(approx.getPoles(True))
+
+ ktar = 4.5
+
+ approx.plotApp(ktar, name = 'u_Pade''')
+ approx.plotHF(ktar, name = 'u_HF')
+ approx.plotErr(ktar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
+ print(appErr, solNorm, appErr/solNorm)
+
+############
+elif testNo == 4:
+
+ def Dboundary(x, on_boundary):
+ return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
+
+ A = 10
+ kappa = 4
+ theta = PI * 90 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ u0ex = - A * sp.exp(1.j * phiex)
+
+ nx = ny = 40
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+
+ DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
+
+ params = {'N':40, 'M':39, 'S':45, 'polyBasis':'CHEBYSHEV', 'POD':True}
+
+ solver = HFEngine.FenicsHelmholtzScatteringAugmentedEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0,
+ FEDegree = 3, DirichletBoundary = Dboundary,
+ RobinBoundary = 'rest', DirichletDatum = DirichletTerm,
+ constraintType = 'IDENTITY')
+ plotter = HSEngine.FenicsHSAugmentedEngine(solver.V, 2)
+ approx = Pade.ROMApproximantLagrangePade(solver, plotter, ks = [0, 8],
+ approxParameters = params,
+ plotSnapshots = False)
+
+ ktar = 4.5
+
+ approx.plotApp(ktar, name = 'u_Pade''')
+ approx.plotHF(ktar, name = 'u_HF')
+ approx.plotErr(ktar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ktar, kappa), approx.HFNorm(ktar, kappa)
+ print(appErr, solNorm, np.divide(appErr, solNorm))
+
+ print(approx.getPoles(True))
diff --git a/examples/Python/HelmholtzPadeTaylorApproximant.py b/examples/Python/HelmholtzPadeTaylorApproximant.py
new file mode 100644
index 0000000..19a52d4
--- /dev/null
+++ b/examples/Python/HelmholtzPadeTaylorApproximant.py
@@ -0,0 +1,196 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Example homogeneous Dirichlet forcing wave
+
+from __future__ import print_function
+import fenics as fen
+import numpy as np
+import sympy as sp
+from context import FenicsHelmholtzEngine as HFEngine
+from context import FenicsHSEngine as HSEngine
+from context import ROMApproximantTaylorPade as Pade
+
+PI = np.pi
+
+testNo = 4
+
+if testNo == 1:
+
+ params = {'N':4, 'M':3, 'E':4, 'sampleType':'Krylov', 'POD':True}
+
+ nu = 12**.5
+ theta = PI / 3
+ z0 = 12+.5j
+ ztar = 11
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
+ phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
+ uex = wex * sp.exp(-1.j * phiex)
+ fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
+
+ nx = ny = 40
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+ forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+ solver = HFEngine.FenicsHelmholtzEngine(mesh = mesh, wavenumber = z0**.5, forcingTerm = forcingTerm, FEDegree = 3,
+ DirichletBoundary = 'all', DirichletDatum = 0)
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = Pade.ROMApproximantTaylorPade(solver, plotter, k0 = z0, w = np.real(z0**.5),
+ approxParameters = params)
+
+ approx.plotApp(ztar, name = 'u_Pade''')
+ approx.plotHF(ztar, name = 'u_HF')
+ approx.plotErr(ztar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
+ print(appErr, solNorm, appErr/solNorm)
+
+ print(approx.getPoles(True))
+
+############
+elif testNo == 2:
+
+ params = {'N':8, 'M':7, 'E':8, 'sampleType':'Arnoldi', 'POD':True}
+ ztar = 3.95 ** 2
+
+ n1 = 2**.5
+ n2 = 3**.5
+ kappa = 4
+ theta = PI * 75 / 180
+ d1, d2 = np.cos(theta), np.sin(theta)
+ K1 = kappa * n1 * d1
+ if kappa * n2 >= K1:
+ K2 = ((kappa*n2)**2 - K1**2)**.5
+ else:
+ K2 = 1.j * (K1**2 - (kappa*n2)**2)**.5
+ R = (kappa * n1 * d2 - K2) / (kappa * n1 * d2 + K2)
+ T = R + 1
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ uex1 = T*sp.exp(1.j*(K1*x+K2*y))
+ uex2 = sp.exp(1.j*kappa*n1*(d1*x+d2*y)) + R*sp.exp(1.j*kappa*n1*(d1*x-d2*y))
+
+ # Exact solution
+ uexRe = fen.Expression('x[1]>=0 ? {0} : {1}'.format(\
+ sp.printing.ccode(sp.re(uex1)), sp.printing.ccode(sp.re(uex2))), degree=4)
+ uexIm = fen.Expression('x[1]>=0 ? {0} : {1}'.format(\
+ sp.printing.ccode(sp.im(uex1)), sp.printing.ccode(sp.im(uex2))), degree=4)
+
+ # wavenumber term
+ nRe = fen.Expression('x[1]<0 ? n1r : n2r',
+ n1r = n1.real, n2r = n2.real, degree=4)
+ nIm = fen.Expression('x[1]<0 ? n1i : n2i',
+ n1i = n1.imag, n2i = n2.imag, degree=4)
+
+ # Create mesh and define function space
+ nx = ny = 50
+ mesh = fen.RectangleMesh(fen.Point(-PI/2,-PI/2), fen.Point(PI/2,PI/2), nx, ny)
+ solver = HFEngine.FenicsHelmholtzEngine(mesh = mesh, wavenumber = kappa, refractionIndex = (nRe, nIm),
+ forcingTerm = 0, FEDegree = 3, DirichletBoundary = 'all',
+ DirichletDatum = (uexRe, uexIm))
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = Pade.ROMApproximantTaylorPade(solver, plotter, k0 = kappa ** 2,
+ w = kappa, approxParameters = params,
+ plotDer = True)
+
+ approx.plotApp(ztar, name = 'u_Pade''')
+ approx.plotHF(ztar, name = 'u_HF')
+ approx.plotErr(ztar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
+ print(appErr, solNorm, appErr/solNorm)
+
+ print(approx.getPoles(True))
+
+############
+elif testNo == 3:
+
+ def Dboundary(x, on_boundary):
+ return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
+
+ A = 10
+ kappa = 2
+ theta = PI * 30 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ wex = 4/PI**2 * x * (PI - x)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ uex = wex * sp.exp(-1.j * phiex)
+ fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
+
+ nx = ny = 50
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+
+ forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+ params = {'N':8, 'M':7, 'E':8, 'sampleType':'Krylov', 'POD':True}
+
+ solver = HFEngine.FenicsHelmholtzScatteringEngine(mesh = mesh, wavenumber = kappa, forcingTerm = forcingTerm,
+ FEDegree = 3, DirichletBoundary = Dboundary,
+ RobinBoundary = 'rest', DirichletDatum = 0)
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = Pade.ROMApproximantTaylorPade(solver, plotter, k0 = kappa,
+ approxParameters = params,
+ plotDer = True)
+
+ approx.setupApprox()
+ print(np.roots(approx.Q[::-1]) + kappa)
+
+ ktar = 1.8 - .3j
+
+ approx.plotApp(ktar, name = 'u_Pade''')
+ approx.plotHF(ktar, name = 'u_HF')
+ approx.plotErr(ktar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
+ print(appErr, solNorm, appErr/solNorm)
+
+ print(approx.getPoles(True))
+
+
+############
+elif testNo == 4:
+
+ def Dboundary(x, on_boundary):
+ return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
+
+ A = 10
+ kappa = 2
+ theta = PI * 30 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ wex = 4/PI**2 * x * (PI - x)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ uex = wex * sp.exp(-1.j * phiex)
+ fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
+
+ nx = ny = 20
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+
+ forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+ params = {'N':8, 'M':7, 'E':8, 'sampleType':'Arnoldi', 'POD':True}
+
+ solver = HFEngine.FenicsHelmholtzScatteringAugmentedEngine(mesh = mesh, wavenumber = kappa, forcingTerm = forcingTerm,
+ FEDegree = 3, DirichletBoundary = Dboundary,
+ RobinBoundary = 'rest', DirichletDatum = 0,
+ constraintType = 'IDENTITY')
+ plotter = HSEngine.FenicsHSAugmentedEngine(solver.V, 2)
+ approx = Pade.ROMApproximantTaylorPade(solver, plotter, k0 = kappa,
+ approxParameters = params,
+ plotDer = True)
+
+ approx.setupApprox()
+ print(np.roots(approx.Q[::-1]) + kappa)
+
+ ktar = 1.8 - .3j
+
+ approx.plotApp(ktar, name = 'u_Pade''')
+ approx.plotHF(ktar, name = 'u_HF')
+ approx.plotErr(ktar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ktar, kappa), approx.HFNorm(ktar, kappa)
+ print(appErr, solNorm, np.divide(appErr, solNorm))
+
+ print(approx.getPoles(True))
diff --git a/examples/Python/HelmholtzRBLagrangeApproximant.py b/examples/Python/HelmholtzRBLagrangeApproximant.py
new file mode 100644
index 0000000..bbb8ebb
--- /dev/null
+++ b/examples/Python/HelmholtzRBLagrangeApproximant.py
@@ -0,0 +1,128 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from __future__ import print_function
+import fenics as fen
+import numpy as np
+import sympy as sp
+from context import FenicsHelmholtzEngine as HFEngine
+from context import FenicsHSEngine as HSEngine
+from context import ROMApproximantLagrangeRB as RB
+
+PI = np.pi
+
+testNo = 3
+
+if testNo == 1:
+
+ params = {'S':5, 'nodesType':'CHEBYSHEV', 'POD':True}
+
+ nu = 12**.5
+ theta = PI / 3
+ ztar = 11
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
+ phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
+ uex = wex * sp.exp(-1.j * phiex)
+ fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
+
+ nx = ny = 40
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+ forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+ solver = HFEngine.FenicsHelmholtzEngine(mesh = mesh, wavenumber = nu, forcingTerm = forcingTerm, FEDegree = 3,
+ DirichletBoundary = 'all', DirichletDatum = 0)
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = RB.ROMApproximantLagrangeRB(solver, plotter, ks = [10 + .5j, 14 + .5j], w = np.real(nu),
+ approxParameters = params)
+
+ approx.plotApp(ztar, name = 'u_RB')
+ approx.plotHF(ztar, name = 'u_HF')
+ approx.plotErr(ztar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
+ print(appErr, solNorm, appErr/solNorm)
+
+ print(approx.getPoles(True))
+
+############
+elif testNo == 2:
+
+ def Dboundary(x, on_boundary):
+ return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
+
+ A = 10
+ kappa = 4
+ theta = PI * 90 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ u0ex = - A * sp.exp(1.j * phiex)
+
+ nx = ny = 40
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+
+ DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
+
+ params = {'S':30, 'nodesType':'CHEBYSHEV', 'POD':True}
+
+ solver = HFEngine.FenicsHelmholtzScatteringEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0,
+ FEDegree = 3, DirichletBoundary = Dboundary,
+ RobinBoundary = 'rest', DirichletDatum = DirichletTerm)
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = RB.ROMApproximantLagrangeRB(solver, plotter, ks = [0, 8],
+ approxParameters = params,
+ plotSnapshots = False)
+
+ ktar = 4.5
+
+ approx.plotApp(ktar, name = 'u_RB')
+ approx.plotHF(ktar, name = 'u_HF')
+ approx.plotErr(ktar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
+ print(appErr, solNorm, appErr/solNorm)
+
+ print(approx.getPoles(True))
+
+############
+elif testNo == 3:
+
+ def Dboundary(x, on_boundary):
+ return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
+
+ A = 10
+ kappa = 4
+ theta = PI * 90 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ u0ex = - A * sp.exp(1.j * phiex)
+
+ nx = ny = 40
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+
+ DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
+
+ params = {'S':30, 'nodesType':'CHEBYSHEV', 'POD':True}
+
+ solver = HFEngine.FenicsHelmholtzScatteringAugmentedEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0,
+ FEDegree = 3, DirichletBoundary = Dboundary,
+ RobinBoundary = 'rest', DirichletDatum = DirichletTerm,
+ constraintType = 'IDENTITY')
+ plotter = HSEngine.FenicsHSAugmentedEngine(solver.V, 2)
+ approx = RB.ROMApproximantLagrangeRB(solver, plotter, ks = [0, 8],
+ approxParameters = params,
+ plotSnapshots = False)
+
+ ktar = 4.5
+
+ approx.plotApp(ktar, name = 'u_RB')
+ approx.plotHF(ktar, name = 'u_HF')
+ approx.plotErr(ktar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ktar, kappa), approx.HFNorm(ktar, kappa)
+ print(appErr, solNorm, np.divide(appErr, solNorm))
+
+ print(approx.getPoles(True))
diff --git a/examples/Python/HelmholtzRBTaylorApproximant.py b/examples/Python/HelmholtzRBTaylorApproximant.py
new file mode 100644
index 0000000..913610a
--- /dev/null
+++ b/examples/Python/HelmholtzRBTaylorApproximant.py
@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Example homogeneous Dirichlet forcing wave
+
+from __future__ import print_function
+import fenics as fen
+import numpy as np
+import sympy as sp
+from context import FenicsHelmholtzEngine as HFEngine
+from context import FenicsHSEngine as HSEngine
+from context import ROMApproximantTaylorRB as RB
+
+PI = np.pi
+
+testNo = 3
+
+if testNo == 1:
+
+ nu = 12**.5
+ theta = PI / 3
+ z0 = 12+.5j
+ ztar = 11
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
+ phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
+ uex = wex * sp.exp(-1.j * phiex)
+ fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
+
+ params = {'E':5, 'POD':True, 'sampleType':'ARNOLDI'}
+
+ nx = ny = 40
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+ forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+ solver = HFEngine.FenicsHelmholtzEngine(mesh = mesh, wavenumber = z0**.5, forcingTerm = forcingTerm, FEDegree = 3,
+ DirichletBoundary = 'all', DirichletDatum = 0)
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = RB.ROMApproximantTaylorRB(solver, plotter, k0 = z0, w = np.real(z0**.5),
+ approxParameters = params)
+
+ approx.plotApp(ztar, name = 'u_RB')
+ approx.plotHF(ztar, name = 'u_HF')
+ approx.plotErr(ztar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
+ print(appErr, solNorm, appErr/solNorm)
+
+ print(approx.getPoles(True))
+
+############
+elif testNo == 2:
+
+ def Dboundary(x, on_boundary):
+ return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
+
+ A = 10
+ kappa = 2
+ theta = PI * 30 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ wex = 4/PI**2 * x * (PI - x)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ uex = wex * sp.exp(-1.j * phiex)
+ fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
+
+ nx = ny = 20
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+
+ forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+ params = {'E':8, 'POD':True}
+
+ solver = HFEngine.FenicsHelmholtzScatteringEngine(mesh = mesh, wavenumber = kappa, forcingTerm = forcingTerm,
+ FEDegree = 3, DirichletBoundary = Dboundary,
+ RobinBoundary = 'rest', DirichletDatum = 0)
+ plotter = HSEngine.FenicsHSEngine(solver.V)
+ approx = RB.ROMApproximantTaylorRB(solver, plotter, k0 = kappa,
+ approxParameters = params,
+ plotDer = False)
+
+ ktar = 1.8 - .2j
+
+ approx.plotApp(ktar, name = 'u_RB')
+ approx.plotHF(ktar, name = 'u_HF')
+ approx.plotErr(ktar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
+ print(appErr, solNorm, appErr/solNorm)
+
+ print(approx.getPoles(True))
+
+############
+elif testNo == 3:
+
+ def Dboundary(x, on_boundary):
+ return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
+
+ A = 10
+ kappa = 2
+ theta = PI * 30 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ wex = 4/PI**2 * x * (PI - x)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ uex = wex * sp.exp(-1.j * phiex)
+ fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
+
+ nx = ny = 20
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+
+ forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+ params = {'E':8, 'POD':True}
+
+ solver = HFEngine.FenicsHelmholtzScatteringAugmentedEngine(mesh = mesh, wavenumber = kappa, forcingTerm = forcingTerm,
+ FEDegree = 3, DirichletBoundary = Dboundary,
+ RobinBoundary = 'rest', DirichletDatum = 0,
+ constraintType = 'MASS')
+ plotter = HSEngine.FenicsHSAugmentedEngine(solver.V, 2)
+ approx = RB.ROMApproximantTaylorRB(solver, plotter, k0 = kappa,
+ approxParameters = params,
+ plotDer = False)
+
+ ktar = 1.8 - .2j
+
+ approx.plotApp(ktar, name = 'u_RB')
+ approx.plotHF(ktar, name = 'u_HF')
+ approx.plotErr(ktar, name = 'err')
+
+ appErr, solNorm = approx.approxError(ktar, kappa), approx.HFNorm(ktar, kappa)
+ print(appErr, solNorm, np.divide(appErr, solNorm))
+
+ print(approx.getPoles(True))
diff --git a/examples/Python/HelmholtzSolver.py b/examples/Python/HelmholtzSolver.py
new file mode 100644
index 0000000..b3386e6
--- /dev/null
+++ b/examples/Python/HelmholtzSolver.py
@@ -0,0 +1,196 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from __future__ import print_function
+import fenics as fen
+import numpy as np
+import sympy as sp
+from context import FenicsHelmholtzEngine as HF
+from context import FenicsHSEngine as HS
+
+testNo = 1
+
+if testNo == 1:
+ PI = np.pi
+
+ def boundary(x, on_boundary):
+ return on_boundary
+
+ nu = 12 ** .5
+ theta = PI / 3
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
+ phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
+ uex = wex * sp.exp(-1.j * phiex)
+ fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
+
+ nx = ny = 40
+ mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI, PI), nx, ny)
+ forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+ solver = HF.FenicsHelmholtzEngine(mesh = mesh, wavenumber = nu + 1.j, forcingTerm = forcingTerm, FEDegree = 3,\
+ DirichletBoundary = boundary, DirichletDatum = 0)
+
+ uh = solver.solve()
+ plotter = HS.FenicsHSEngine(solver.V)
+ print(plotter.norm(uh, np.real(nu)))
+ plotter.plot(uh)
+
+
+###########
+elif testNo == 2:
+
+ def boundary(x, on_boundary):
+ return on_boundary
+
+ PI = np.pi
+ n1 = 4**.5
+ n2 = 1**.5
+ kappa = 3
+ theta = PI * 70 / 180
+ d1, d2 = np.cos(theta), np.sin(theta)
+ K1 = kappa * n1 * d1
+ if kappa * n2 >= K1:
+ K2 = ((kappa*n2)**2 - K1**2)**.5
+ else:
+ K2 = 1.j * (K1**2 - (kappa*n2)**2)**.5
+ R = (kappa * n1 * d2 - K2) / (kappa * n1 * d2 + K2)
+ T = R + 1
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ uex1 = T*sp.exp(1.j*(K1*x+K2*y))
+ uex2 = sp.exp(1.j*kappa*n1*(d1*x+d2*y)) + R*sp.exp(1.j*kappa*n1*(d1*x-d2*y))
+
+ # Exact solution
+ uexRe = fen.Expression('x[1]>=0 ? {0} : {1}'.format(\
+ sp.printing.ccode(sp.re(uex1)), sp.printing.ccode(sp.re(uex2))), degree=4)
+ uexIm = fen.Expression('x[1]>=0 ? {0} : {1}'.format(\
+ sp.printing.ccode(sp.im(uex1)), sp.printing.ccode(sp.im(uex2))), degree=4)
+
+ # refraction index
+ n2Re = fen.Expression('x[1]<0 ? n1r : n2r',
+ n1r = n1.real, n2r = n2.real, degree=4)
+ n2Im = fen.Expression('x[1]<0 ? n1i : n2i',
+ n1i = n1.imag, n2i = n2.imag, degree=4)
+
+ # Create mesh and define function space
+ nx = ny = 50
+ mesh = fen.RectangleMesh(fen.Point(-PI/2,-PI/2), fen.Point(PI/2,PI/2), nx, ny)
+
+ solver = HF.FenicsHelmholtzEngine(mesh = mesh, wavenumber = kappa, refractionIndex = (n2Re, n2Im),\
+ forcingTerm = 0, FEDegree = 3, DirichletBoundary = boundary,\
+ DirichletDatum = (uexRe, uexIm))
+
+ uh = solver.solve()
+ plotter = HS.FenicsHSEngine(solver.V)
+ print(plotter.norm(uh, kappa))
+ plotter.plot(uh)
+
+
+###########
+elif testNo == 3:
+
+ import mshr
+ from matplotlib import pyplot as plt
+
+ PI = np.pi
+ R = 5
+ def Dboundary(x, on_boundary):
+ return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
+
+ A = 10
+ kappa = 12**.5
+ theta = - PI * 90 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ u0ex = - A * sp.exp(1.j * phiex)
+
+ npoints = 50
+ scatterer = mshr.Polygon([fen.Point(-1, -.5),
+ fen.Point(1, -.5),
+ fen.Point(1, .5),
+ fen.Point(.8, .5),
+ fen.Point(.8, -.3),
+ fen.Point(-.8, -.3),
+ fen.Point(-.8, .5),
+ fen.Point(-1, .5),
+ ])
+ mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0, 0), R) - scatterer, npoints)
+
+ plt.jet()
+ plt.figure()
+ fen.plot(mesh)
+
+ DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
+
+ solver = HF.FenicsHelmholtzScatteringEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0, FEDegree = 3,\
+ DirichletBoundary = Dboundary, RobinBoundary = 'rest',\
+ DirichletDatum = DirichletTerm)
+
+ baseRe, baseIm = DirichletTerm
+ baseRe = fen.project(fen.Expression(baseRe, degree = 4), solver.V)
+ baseIm = fen.project(fen.Expression(baseIm, degree = 4), solver.V)
+ uinc = np.array(baseRe.vector()) + 1.j * np.array(baseIm.vector())
+
+ uh = solver.solve()
+ plotter = HS.FenicsHSEngine(solver.V)
+ print(plotter.norm(uh, kappa))
+ print(plotter.norm(uh - uinc, kappa))
+ plotter.plot(uh)
+ plotter.plot(uh - uinc)
+
+###########
+elif testNo == 4:
+
+ import mshr
+ from matplotlib import pyplot as plt
+
+ PI = np.pi
+ R = 5
+ def Dboundary(x, on_boundary):
+ return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
+
+ A = 10
+ kappa = 12**.5
+ theta = - PI * 90 / 180
+
+ x, y = sp.symbols('x[0] x[1]', real=True)
+ phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
+ u0ex = - A * sp.exp(1.j * phiex)
+
+ npoints = 40
+ scatterer = mshr.Polygon([fen.Point(-1, -.5),
+ fen.Point(1, -.5),
+ fen.Point(1, .5),
+ fen.Point(.8, .5),
+ fen.Point(.8, -.3),
+ fen.Point(-.8, -.3),
+ fen.Point(-.8, .5),
+ fen.Point(-1, .5),
+ ])
+ mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0, 0), R) - scatterer, npoints)
+
+ plt.jet()
+ plt.figure()
+ fen.plot(mesh)
+
+ DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))), sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
+
+ solver = HF.FenicsHelmholtzScatteringAugmentedEngine(mesh = mesh, wavenumber = kappa, forcingTerm = 0, FEDegree = 3,\
+ DirichletBoundary = Dboundary, RobinBoundary = 'rest',\
+ DirichletDatum = DirichletTerm, constraintType = 'MASS')
+
+ baseRe, baseIm = DirichletTerm
+ baseRe = fen.project(fen.Expression(baseRe, degree = 4), solver.V)
+ baseIm = fen.project(fen.Expression(baseIm, degree = 4), solver.V)
+ uinc = np.array(baseRe.vector()) + 1.j * np.array(baseIm.vector())
+ uinc = np.concatenate((uinc, kappa * uinc))
+
+ uh = solver.solve()
+ plotter = HS.FenicsHSAugmentedEngine(solver.V, 2)
+ print(plotter.norm(uh, kappa))
+ print(plotter.norm(uh - uinc, kappa))
+ plotter.plot(uh)
+ plotter.plot(uh - uinc)
diff --git a/examples/Python/HelmholtzTaylorApproximantsSweep.py b/examples/Python/HelmholtzTaylorApproximantsSweep.py
new file mode 100644
index 0000000..c4edf75
--- /dev/null
+++ b/examples/Python/HelmholtzTaylorApproximantsSweep.py
@@ -0,0 +1,100 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Example homogeneous Dirichlet forcing wave SWEEP
+
+from __future__ import print_function
+import fenics as fen
+import numpy as np
+import sympy as sp
+from context import FenicsHelmholtzEngine as HFEngine
+from context import FenicsHSEngine as HSEngine
+from context import ROMApproximantTaylorPade as Pade
+from context import ROMApproximantTaylorRB as RB
+from context import ROMApproximantSweeper as Sweeper
+
+PI = np.pi
+
+nu = 12**.5
+theta = PI / 3
+z0 = 12 + .5j
+npoints = 31
+ktars = np.linspace(7, 16, npoints)
+
+x, y = sp.symbols('x[0] x[1]', real=True)
+wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
+phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
+uex = wex * sp.exp(-1.j * phiex)
+fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
+
+nx = ny = 10
+mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+solver = HFEngine.FenicsHelmholtzEngine(mesh = mesh, wavenumber = nu, forcingTerm = forcingTerm, FEDegree = 3,
+ DirichletBoundary = 'all', DirichletDatum = 0)
+plotter = HSEngine.FenicsHSEngine(solver.V)
+
+shift = 5
+nsets = 5
+stride = 2
+Emax = stride * (nsets - 1) + shift + 2
+params = {'Emax':Emax, 'sampleType':'ARNOLDI', 'POD':True}
+paramsSetsPade = [None] * nsets
+paramsSetsRB = [None] * nsets
+for i in range(nsets):
+ paramsSetsPade[i] = {'N':stride*i+shift+1, 'M':stride*i+shift,
+ 'E':stride*i+shift+1}
+ paramsSetsRB[i] = {'E':stride*i+shift+1,'R':stride*i+shift+2}
+
+appPade = Pade.ROMApproximantTaylorPade(solver, plotter, k0 = z0, w = np.real(z0**.5),
+ approxParameters = params)
+appRB = RB.ROMApproximantTaylorRB(solver, plotter, k0 = z0, w = np.real(z0**.5),
+ approxParameters = params)
+
+sweeper = Sweeper.ROMApproximantSweeper(ktars = ktars, mostExpensive = 'Approx')
+sweeper.ROMEngine = appPade
+sweeper.params = paramsSetsPade
+filenamePade = sweeper.sweep('../Data/HelmholtzBubbleTaylorPade.dat', outputs = 'ALL')
+
+sweeper.ROMEngine = appRB
+sweeper.params = paramsSetsRB
+filenameRB = sweeper.sweep('../Data/HelmholtzBubbleTaylorRB.dat', outputs = 'ALL')
+
+####################
+
+from matplotlib import pyplot as plt
+plt.jet()
+
+for i in range(nsets):
+ nDerivatives = stride*i+shift+1
+ PadeOutput = sweeper.read(filenamePade, {'E':nDerivatives},
+ ['kRe', 'HFNorm', 'AppNorm', 'AppError'])
+ RBOutput = sweeper.read(filenameRB, {'E':nDerivatives},
+ ['kRe', 'AppNorm', 'AppError'])
+
+ ktarsF = PadeOutput['kRe']
+ solNormF = PadeOutput['HFNorm']
+ PadektarsF = PadeOutput['kRe']
+ PadeNormF = PadeOutput['AppNorm']
+ PadeErrorF = PadeOutput['AppError']
+ RBktarsF = RBOutput['kRe']
+ RBNormF = RBOutput['AppNorm']
+ RBErrorF = RBOutput['AppError']
+
+ plt.figure()
+ plt.semilogy(ktarsF, solNormF, 'k-', label='Sol norm')
+ plt.semilogy(PadektarsF, PadeNormF, 'b.--', label='Pade'' norm, E = {}'.format(nDerivatives))
+ plt.semilogy(RBktarsF, RBNormF, 'g.--', label='RB norm, E = {}'.format(nDerivatives))
+ plt.legend()
+ plt.grid()
+ plt.figure()
+ plt.semilogy(PadektarsF, PadeErrorF, 'b', label='Pade'' error, E = {}'.format(nDerivatives))
+ plt.semilogy(RBktarsF, RBErrorF, 'g', label='RB error, E = {}'.format(nDerivatives))
+ plt.legend()
+ plt.grid()
+ plt.figure()
+ plt.semilogy(ktarsF, PadeErrorF / solNormF, 'b', label='Pade'' relative error, E = {}'.format(nDerivatives))
+ plt.semilogy(RBktarsF, RBErrorF / solNormF, 'g', label='RB relative error, E = {}'.format(nDerivatives))
+ plt.legend()
+ plt.grid()
diff --git a/examples/Python/HelmholtzTaylorPoleIdentification.py b/examples/Python/HelmholtzTaylorPoleIdentification.py
new file mode 100644
index 0000000..c302515
--- /dev/null
+++ b/examples/Python/HelmholtzTaylorPoleIdentification.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Example homogeneous Dirichlet forcing wave
+
+from __future__ import print_function
+import fenics as fen
+import numpy as np
+import sympy as sp
+from context import utilities
+from context import FenicsHelmholtzEngine as HFEngine
+from context import FenicsHSEngine as HSEngine
+from context import ROMApproximantTaylorPade as Pade
+from context import ROMApproximantTaylorRB as RB
+
+PI = np.pi
+
+nu = 12**.5
+z0 = 12+0.j
+theta = PI / 3
+
+x, y = sp.symbols('x[0] x[1]', real=True)
+wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
+phiex = nu * (x * np.cos(theta) + y * np.sin(theta))
+uex = wex * sp.exp(-1.j * phiex)
+fex = - uex.diff(x, 2) - uex.diff(y, 2) - nu**2 * uex
+
+nx = ny = 25
+mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI,PI), nx, ny)
+forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))), sp.printing.ccode(sp.simplify(sp.im(fex)))]
+
+Nmin, Nmax = 2, 10
+Nvals = np.arange(Nmin, Nmax + 1, 2)
+
+params = {'N':Nmin, 'M':0, 'Emax':Nmax, 'POD':True, 'sampleType':'Arnoldi'}#, 'robustTol':1e-14}
+
+#boolCon = lambda x : np.abs(np.imag(x)) < 1e-1 * np.abs(np.real(x) - np.real(z0))
+#cleanupParameters = {'boolCondition':boolCon, 'residueCheck':True}
+
+solver = HFEngine.FenicsHelmholtzEngine(mesh = mesh, wavenumber = z0**.5, forcingTerm = forcingTerm, FEDegree = 3,
+ DirichletBoundary = 'all', DirichletDatum = 0)
+plotter = HSEngine.FenicsHSEngine(solver.V)
+approxP = Pade.ROMApproximantTaylorPade(solver, plotter, k0 = z0, w = np.real(z0**.5),
+ approxParameters = params)#, equilibration = True,
+# cleanupParameters = cleanupParameters)
+approxR = RB.ROMApproximantTaylorRB(solver, plotter, k0 = z0, w = np.real(z0**.5),
+ approxParameters = params)
+
+rP, rE = [None] * len(Nvals), [None] * len(Nvals)
+
+verbose = 1
+for j, N in enumerate(Nvals):
+ if verbose > 0:
+ print('N = E = {}'.format(N))
+ approxP.approxParameters = {'N':N, 'E':N}
+ approxR.approxParameters = {'R':N, 'E':N}
+ if verbose > 1:
+ print(approxP.approxParameters)
+ print(approxR.approxParameters)
+
+ rP[j] = approxP.getPoles(True)
+ rE[j] = approxR.getPoles(True)
+ if verbose > 2:
+ print(rP)
+ print(rE)
+
+from matplotlib import pyplot as plt
+plt.set_cmap('jet')
+plotRows = int(np.ceil(len(Nvals) / 3))
+fig, axes = plt.subplots(plotRows, 3, figsize = (15, 3.5 * plotRows))
+for j, N in enumerate(Nvals):
+ i1, i2 = int(np.floor(j / 3)), j % 3
+ axes[i1, i2].set_title('N = E = {}'.format(N))
+ axes[i1, i2].plot(np.real(rP[j]), np.imag(rP[j]), 'Xb',
+ label="Pade'", markersize = 8)
+ axes[i1, i2].plot(np.real(rE[j]), np.imag(rE[j]), '*r',
+ label="RB", markersize = 10)
+ axes[i1, i2].axhline(linewidth=1, color='k')
+ xmin, xmax = axes[i1, i2].get_xlim()
+ res = utilities.squareResonances(xmin, xmax, False)
+ axes[i1, i2].plot(res, np.zeros_like(res), 'ok', markersize = 4)
+ axes[i1, i2].grid()
+ axes[i1, i2].set_xlim(xmin, xmax)
+ axes[i1, i2].axis('equal')
+ p = axes[i1, i2].legend()
+plt.tight_layout()
+for j in range((len(Nvals) - 1) % 3 + 1, 3):
+ axes[plotRows - 1, j].axis('off')
diff --git a/examples/Python/__init__.py b/examples/Python/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/examples/Python/context.py b/examples/Python/context.py
new file mode 100644
index 0000000..70ddf70
--- /dev/null
+++ b/examples/Python/context.py
@@ -0,0 +1,11 @@
+# -*- coding: utf-8 -*-
+
+import sys, os
+module_path = os.path.abspath(os.path.join(os.path.dirname(__file__), '..', '..', 'main'))
+if module_path not in sys.path:
+ sys.path.append(module_path)
+
+import utilities, ROMApproximant, ROMApproximantSweeper
+import FenicsHSEngine, FenicsHelmholtzEngine, FreeFemHSEngine, FreeFemHelmholtzEngine, FreeFemConversionTools
+import ROMApproximantTaylor, ROMApproximantTaylorPade, ROMApproximantTaylorRB
+import ROMApproximantLagrange, ROMApproximantLagrangePade, ROMApproximantLagrangeRB
diff --git a/examples/Python/test.py b/examples/Python/test.py
new file mode 100644
index 0000000..577aeed
--- /dev/null
+++ b/examples/Python/test.py
@@ -0,0 +1,17 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from __future__ import print_function
+import numpy as np
+from context import FreeFemHSEngine as FFHSE
+
+V = ("load \"Element_P3\";\n"
+ "int n = 5;\n"
+ "mesh Th = square(n, n, [pi * x, pi * y]);\n"
+ "fespace V(Th, P3);")
+
+u = np.random.rand(256)
+
+engine = FFHSE.FreeFemHSEngine(V, 2, "Th", "V")
+
+nm = engine.norm(u, 'H1')
\ No newline at end of file
diff --git a/examples/__init__.py b/examples/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/examples/context.py b/examples/context.py
new file mode 100644
index 0000000..2fdc3c4
--- /dev/null
+++ b/examples/context.py
@@ -0,0 +1,11 @@
+# -*- coding: utf-8 -*-
+
+import sys, os
+module_path = os.path.abspath(os.path.join(os.path.dirname(__file__), '..', 'main'))
+if module_path not in sys.path:
+ sys.path.append(module_path)
+
+import utilities, ROMApproximant, ROMApproximantSweeper
+import FenicsHSEngine, FenicsHelmholtzEngine, FreeFemHSEngine, FreeFemHelmholtzEngine, FreeFemConversionTools
+import ROMApproximantTaylor, ROMApproximantTaylorPade, ROMApproximantTaylorRB
+import ROMApproximantLagrange, ROMApproximantLagrangePade, ROMApproximantLagrangeRB
diff --git a/main/FenicsHSEngine.py b/main/FenicsHSEngine.py
new file mode 100644
index 0000000..301da21
--- /dev/null
+++ b/main/FenicsHSEngine.py
@@ -0,0 +1,202 @@
+#!/usr/bin/python
+
+import numpy as np
+from matplotlib import pyplot as plt
+import fenics as fen
+
+class FenicsHSEngine:
+ """
+ Fenics-based Hilbert space engine.
+ """
+
+ def __init__(self, V:"Fenics FE space"):
+ self.V = V
+
+ def name(self) -> str:
+ """Class label."""
+ return self.__class__.__name__
+
+ def norm(self, u:"numpy 1D array",
+ normType : "numpy 2D array, number or str" = "H10") -> float:
+ """
+ Compute general norm of complex-valued function with given dofs.
+
+ Args:
+ u: numpy complex array with function dofs.
+ normType(optional): Target norm identifier. If matrix, target norm
+ is one induced by normType. If number, target norm is weighted
+ H^1 norm with given weight. If string, must be recognizable by
+ Fenics norm command. Defaults to 'H10'.
+
+ Returns:
+ Norm of the function (non-negative).
+ """
+ if type(normType).__name__[-6:] == "matrix":
+ return np.abs(u.dot(normType.dot(u).conj())) ** .5
+ if isinstance(normType, (int, float)):
+ return (FenicsHSEngine.norm(self, u, "H10")**2
+ + (normType * FenicsHSEngine.norm(self, u, "L2"))**2)**.5
+ uRe = fen.Function(self.V)
+ uIm = fen.Function(self.V)
+ uRe.vector()[:] = np.array(np.real(u), dtype = float)
+ uIm.vector()[:] = np.array(np.imag(u), dtype = float)
+ return (fen.norm(uRe, normType)**2 + fen.norm(uIm, normType)**2)**.5
+
+ def plot(self, u:"numpy 1D array", name : str = "u", **figspecs):
+ """
+ Do some nice plots of the complex-valued function with given dofs.
+
+ Args:
+ u: numpy complex array with function dofs.
+ name(optional): Name to be shown as title of the plots. Defaults to
+ 'u'.
+ figspecs(optional key args): Optional arguments for matplotlib
+ figure creation.
+ """
+ if 'figsize' not in figspecs.keys(): figspecs['figsize'] = (13, 3)
+
+ uRe = fen.Function(self.V)
+ uIm = fen.Function(self.V)
+ uRe.vector()[:] = np.array(np.real(u), dtype = float)
+ uIm.vector()[:] = np.array(np.imag(u), dtype = float)
+
+ uAbs = fen.project(np.power(np.power(uRe, 2) + np.power(uIm, 2),
+ .5), self.V)
+ uPhase = fen.project(fen.atan(np.divide(uIm, uRe)), self.V)
+ combined_data_ReIm = np.concatenate([np.real(u), np.imag(u)])
+ _min, _max = np.amin(combined_data_ReIm), np.amax(combined_data_ReIm)
+
+ plt.figure(**figspecs)
+
+ plt.jet()
+
+ plt.subplot(141)
+ p = fen.plot(uAbs, title = "|{0}|".format(name))
+ plt.colorbar(p)
+
+ plt.subplot(142)
+ p = fen.plot(uPhase, title = "phase({0})".format(name))
+ plt.colorbar(p)
+
+ plt.subplot(143)
+ p = fen.plot(uRe, title = "Re({0})".format(name),
+ vmin = _min, vmax = _max)
+ plt.colorbar(p)
+
+ plt.subplot(144)
+ p = fen.plot(uIm, title = "Im({0})".format(name),
+ vmin = _min, vmax = _max)
+ plt.colorbar(p)
+ plt.tight_layout()
+ plt.show()
+
+
+class FenicsHSAugmentedEngine(FenicsHSEngine):
+ """
+ Fenics-based Hilbert space engine for augmented spaces.
+ """
+ def __init__(self, V:"Fenics FE space", d : int = 1):
+ self.d = d
+ FenicsHSEngine.__init__(self, V)
+
+ def getActualSize(self, u:"numpy 1D array"):
+ """
+ Compute size of unaugmented vector batches.
+
+ Args:
+ u: numpy complex array with function dofs.
+
+ Returns:
+ Batch size.
+ """
+ N = int(len(u) / self.d)
+ if not np.isclose(len(u), self.d * N):
+ raise Exception(("Input vector lenght must be multiple of"
+ "augmentation dimension d."))
+ return N
+
+ def norm(self, u:"numpy 1D array",
+ normType : "numpy 2D array, number or str" = "H10")\
+ -> "numpy 1D array":
+ """
+ Compute general norm of complex-valued function with given dofs.
+
+ Args:
+ u: numpy complex array with function dofs.
+ normType(optional): Target norm identifier. If matrix, target norm
+ is one induced by normType. If number, target norm is weighted
+ H^1 norm with given weight. If string, must be recognizable by
+ Fenics norm command. Defaults to 'H10'.
+
+ Returns:
+ Norms of the function (non-negative).
+ """
+ N = self.getActualSize(u)
+ if isinstance(normType, (int, float)):
+ norms = [None] * self.d
+ for j in range(self.d):
+ uj = u[j * N : (j + 1) * N]
+ norms[j] = FenicsHSEngine.norm(self, uj, normType)
+ return norms
+ else:
+ if (type(normType).__name__[-6:] == "matrix"
+ and normType.shape[0] % N == 0):
+ N = normType.shape[0]
+ else:
+ raise Exception("Energy matrix dimension mismatch.")
+ return FenicsHSEngine.norm(self, u[: N], normType)
+
+ def plot(self, u:"numpy 1D array", name : str = "u", **figspecs):
+ """
+ Do some nice plots of the complex-valued function with given dofs.
+
+ Args:
+ u: numpy complex array with function dofs.
+ name(optional): Name to be shown as title of the plots. Defaults to
+ 'u'.
+ figspecs(optional key args): Optional arguments for matplotlib
+ figure creation.
+ """
+ if 'figsize' not in figspecs.keys(): figspecs['figsize'] = (13, 3)
+
+ N = self.getActualSize(u)
+
+ for j in range(self.d):
+ uj = u[j * N : (j + 1) * N]
+
+ uRe = fen.Function(self.V)
+ uIm = fen.Function(self.V)
+ uRe.vector()[:] = np.array(np.real(uj), dtype = float)
+ uIm.vector()[:] = np.array(np.imag(uj), dtype = float)
+
+ uAbs = fen.project(np.power(np.power(uRe, 2) + np.power(uIm, 2),
+ .5), self.V)
+ uPhase = fen.project(fen.atan(np.divide(uIm, uRe)), self.V)
+ combined_data_ReIm = np.concatenate([np.real(uj), np.imag(uj)])
+ _min = np.amin(combined_data_ReIm)
+ _max = np.amax(combined_data_ReIm)
+
+ plt.figure(**figspecs)
+
+ plt.jet()
+
+ plt.subplot(141)
+ p = fen.plot(uAbs, title = "|{0}|, comp. {1}".format(name, j))
+ plt.colorbar(p)
+
+ plt.subplot(142)
+ p = fen.plot(uPhase,
+ title = "phase({0}, comp. {1})".format(name, j))
+ plt.colorbar(p)
+
+ plt.subplot(143)
+ p = fen.plot(uRe, title = "Re({0}, comp. {1})".format(name, j),
+ vmin = _min, vmax = _max)
+ plt.colorbar(p)
+
+ plt.subplot(144)
+ p = fen.plot(uIm, title = "Im({0}, comp. {1})".format(name, j),
+ vmin = _min, vmax = _max)
+ plt.colorbar(p)
+ plt.tight_layout()
+ plt.show()
diff --git a/main/FenicsHelmholtzEngine.py b/main/FenicsHelmholtzEngine.py
new file mode 100644
index 0000000..255e4f5
--- /dev/null
+++ b/main/FenicsHelmholtzEngine.py
@@ -0,0 +1,958 @@
+#!/usr/bin/python
+
+from copy import copy
+import warnings
+import numpy as np
+import sympy as sp
+import sympy.printing as sprint
+import scipy.sparse as scsp
+import scipy.sparse.linalg as spla
+import fenics as fen
+
+PI = np.pi
+fenZERO = fen.Constant(0)
+fenZEROC = [fenZERO, fenZERO]
+
+class DomainError(Exception):
+ """
+ Domain error exception.
+
+ Args:
+ value: Human readable string describing the exception.
+
+ Attributes:
+ value: Human readable string describing the exception.
+ """
+ def __init__(self, value:str):
+ self.value = value
+ def __str__(self):
+ return self.value
+
+def CustomExpressionParser(value:"expression",
+ degree:int) -> "Fenics Expression":
+ """
+ Numerical and Fenics expressions parser.
+
+ Args:
+ value: Expression to be parsed. Accepts 2-tuples composed of real and
+ imaginary parts. Available elementary formats are numbers, strings
+ and Fenics Expression's.
+ degree: Degree of Fenics FE interpolant of expression.
+
+ Returns:
+ Fenics FE interpolant of expression.
+ """
+ try:
+ if isinstance(value, (list, tuple)):
+ if len(value) == 1:
+ return CustomExpressionParser(value[0])
+ elif len(value) != 2:
+ raise Exception("Parsing error")
+ if isinstance(value[0], (int, float, complex)):
+ if any([isinstance(y, complex) for y in value]):
+ raise Exception("Parsing error")
+ valueRe = fen.Constant(value[0])
+ valueIm = fen.Constant(value[1])
+ elif isinstance(value[0], str):
+ x, y, z, x0, x1, x2 = sp.symbols("x y z x[0] x[1] x[2]",
+ real=True)
+ xyDict = {"x": x, "y": y, "z": z}
+ valueReParsed = value[0].replace("x[0]", "x")\
+ .replace("x[1]", "y")\
+ .replace("x[2]", "z")
+ valueImParsed = value[1].replace("x[0]", "x")\
+ .replace("x[1]", "y")\
+ .replace("x[2]", "z")
+ valueReSym = sp.sympify(valueReParsed, locals=xyDict)\
+ .subs([(x, x0), (y, x1), (z, x2)])
+ valueImSym = sp.sympify(valueImParsed, locals=xyDict)\
+ .subs([(x, x0), (y, x1), (z, x2)])
+ valueReStr = sprint.ccode(valueReSym).rpartition("\n")[-1]
+ valueImStr = sprint.ccode(valueImSym).rpartition("\n")[-1]
+ valueRe = fen.Expression(valueReStr, degree = degree)
+ valueIm = fen.Expression(valueImStr, degree = degree)
+ else:
+ valueRe = value[0]
+ valueIm = value[1]
+ else:
+ if isinstance(value, (int, float, complex)):
+ valueRe = fen.Constant(np.real(value))
+ valueIm = fen.Constant(np.imag(value))
+ elif isinstance(value, str):
+ x, y, z, x0, x1, x2 = sp.symbols("x y z x[0] x[1] x[2]",
+ real=True)
+ xyDict = {"x": x, "y": y, "z": z}
+ valueParsed = value.replace("x[0]", "x").replace("x[1]", "y")\
+ .replace("x[2]", "z")
+ valueSym = sp.sympify(valueParsed, locals=xyDict)\
+ .subs([(x, x0), (y, x1), (z, x2)])
+ valueReStr = sprint.ccode(sp.re(valueSym)).rpartition("\n")[-1]
+ valueImStr = sprint.ccode(sp.im(valueSym)).rpartition("\n")[-1]
+ valueImStr = sprint.ccode(valueImSym).rpartition("\n")[-1]
+ valueRe = fen.Expression(valueReStr, degree = degree)
+ valueIm = fen.Expression(valueImStr, degree = degree)
+ else:
+ valueRe = value
+ valueIm = fenZERO
+ except:
+ raise Exception("Parsing error")
+ return copy(valueRe), copy(valueIm)
+
+class FenicsHelmholtzEngine:
+ """
+ Fenics-based solver for Helmholtz problems.
+ - \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega
+ u = u0 on \Gamma_D
+ \partial_nu = g1 on \Gamma_N
+ \partial_nu + h u = g2 on \Gamma_R
+
+ Args:
+ mesh: Domain of Helmholtz problem.
+ FEDegree: FE degree.
+ wavenumber2: Value of k^2.
+ diffusivity(optional): Value of a. Defaults to identically 1.
+ refractionIndex(optional): Value of n. Defaults to identically 1.
+ forcingTerm(optional): Value of f. Defaults to identically 0.
+ DirichletBoundary(optional): Function flagging Dirichlet boundary as
+ True, in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are
+ accepted. Defaults to False everywhere.
+ NeumannBoundary(optional): Function flagging Neumann boundary as True,
+ in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
+ Defaults to False everywhere.
+ RobinBoundary(optional): Function flagging Robin boundary as True, in
+ Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
+ Defaults to False everywhere.
+ DirichletDatum(optional): Value of u0. Defaults to identically 0.
+ NeumannDatum(optional): Value of g1. Defaults to identically 0.
+ RobinDatum(optional): Value of h. Defaults to identically 0.
+ RobinDatum2(optional): Value of g2. Defaults to identically 0.
+
+ Attributes:
+ mesh: Domain of Helmholtz problem.
+ FEDegree: FE degree.
+ wavenumber2: Copy of processed wavenumber2 parameter.
+ diffusivity: Copy of processed diffusivity parameter.
+ refractionIndex: Copy of processed refractionIndex parameter.
+ forcingTerm: Copy of processed forcingTerm parameter.
+ DirichletBoundary: Copy of processed DirichletBoundary parameter.
+ NeumannBoundary: Copy of processed NeumannBoundary parameter.
+ RobinBoundary: Copy of processed RobinBoundary parameter.
+ DirichletDatum: Copy of processed DirichletDatum parameter.
+ NeumannDatum: Copy of processed NeumannDatum parameter.
+ RobinDatum: Copy of processed RobinDatum parameter.
+ RobinDatum2: Copy of processed RobinDatum2 parameter.
+ LU: Whether to use LU solver for computation of derivatives.
+ V: Real FE space.
+ u: Helmholtz problem solution as numpy complex vector.
+ v1: Generic trial functions for variational form evaluation.
+ v2: Generic test functions for variational form evaluation.
+ ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
+ nRe,nIm: Real and imaginary parts of n.
+ n2Re,n2Im: Real and imaginary parts of n^2.
+ fRe,fIm: Real and imaginary parts of f.
+ u0Re,u0Im: Real and imaginary parts of u0.
+ g1Re,g1Im: Real and imaginary parts of g1.
+ hRe,hIm: Real and imaginary parts of h.
+ g2Re,g2Im: Real and imaginary parts of g2.
+ DirichletBCRe,DirichletBCIm: Fenics BC manager for real and imaginary
+ parts of Dirichlet data.
+ A*: Scipy sparse array representation (in CSC format) of A*.
+ b*Re,b*Im: Numpy array representation of real and imaginary parts of
+ b*.
+ invA: Numpy LU solver for A.
+ """
+
+ def __init__(self, mesh:"Fenics mesh",
+ FEDegree:int,
+ wavenumber:"custom complex",
+ diffusivity : "custom expression" = 1,
+ refractionIndex : "custom expression" = 1,
+ forcingTerm : "custom expression" = fenZEROC,
+ DirichletBoundary : "bool function or str" = None,
+ NeumannBoundary : "bool function or str" = None,
+ RobinBoundary : "bool function or str" = None,
+ DirichletDatum : "custom expression" = fenZEROC,
+ NeumannDatum : "custom expression" = fenZEROC,
+ RobinDatum : "custom expression" = fenZEROC,
+ RobinDatum2 : "custom expression" = fenZEROC):
+ self.mesh = copy(mesh)
+ self.FEDegree = FEDegree
+
+ # process boundaries
+ boundariesList = [DirichletBoundary, NeumannBoundary, RobinBoundary]
+ for i in range(3):
+ if isinstance(boundariesList[i], str):
+ boundariesList[i] = boundariesList[i].upper()
+ if boundariesList[i] == "NONE": boundariesList[i] = None
+ DirichletBoundary, NeumannBoundary, RobinBoundary = boundariesList
+ if boundariesList.count(None) == 3:
+ raise DomainError("At least one boundary must be prescribed.")
+ if boundariesList.count("ALL") + boundariesList.count("REST") >= 2:
+ raise DomainError("Only one keyword 'ALL'/'REST' can be used.")
+ def DirichletB(x, on_boundary):
+ if DirichletBoundary is None:
+ return False
+ elif DirichletBoundary == "ALL":
+ return on_boundary
+ elif DirichletBoundary == "REST":
+ return (on_boundary
+ and not self.NeumannBoundary(x, on_boundary)
+ and not self.RobinBoundary(x, on_boundary))
+ else:
+ return DirichletBoundary(x, on_boundary)
+ def NeumannB(x, on_boundary):
+ if NeumannBoundary is None:
+ return False
+ elif NeumannBoundary == "ALL":
+ return on_boundary
+ elif NeumannBoundary == "REST":
+ return (on_boundary
+ and not self.DirichletBoundary(x, on_boundary)
+ and not self.RobinBoundary(x, on_boundary))
+ else:
+ return NeumannBoundary(x, on_boundary)
+ def RobinB(x, on_boundary):
+ if RobinBoundary is None:
+ return False
+ elif RobinBoundary == "ALL":
+ return on_boundary
+ elif RobinBoundary == "REST":
+ return (on_boundary
+ and not self.DirichletBoundary(x, on_boundary)
+ and not self.NeumannBoundary(x, on_boundary))
+ else:
+ return RobinBoundary(x, on_boundary)
+ self.DirichletBoundary = copy(DirichletB)
+ self.NeumannBoundary = copy(NeumannB)
+ self.RobinBoundary = copy(RobinB)
+
+ # create boundary measure
+ class NBoundary(fen.SubDomain):
+ def inside(self, x, on_boundary):
+ return NeumannB(x, on_boundary)
+ class RBoundary(fen.SubDomain):
+ def inside(self, x, on_boundary):
+ return RobinB(x, on_boundary)
+ boundary_markers = fen.FacetFunction("size_t", self.mesh)
+ NBoundary().mark(boundary_markers, 0)
+ RBoundary().mark(boundary_markers, 1)
+ self.ds = fen.Measure("ds", domain=mesh,
+ subdomain_data=boundary_markers)
+
+ self.V = fen.FunctionSpace(self.mesh, "P", self.FEDegree)
+ self.v1 = fen.TrialFunction(self.V)
+ self.v2 = fen.TestFunction(self.V)
+ self.wavenumber2 = wavenumber ** 2
+ self.diffusivity = diffusivity
+ self.refractionIndex = refractionIndex
+ self.forcingTerm = forcingTerm
+ self.DirichletDatum = DirichletDatum
+ self.NeumannDatum = NeumannDatum
+ self.RobinDatum = RobinDatum
+ self.RobinDatum2 = RobinDatum2
+
+ def energyNormMatrix(self, w:float) -> "CSC sparse matrix":
+ """
+ Sparse matrix (in CSR format) assoociated to w-weighted H10 norm.
+
+ Args:
+ w: Weight.
+
+ Returns:
+ Sparse matrix in CSR format.
+ """
+ normMatFen = fen.assemble((
+ fen.dot(fen.grad(self.v1), fen.grad(self.v2))
+ + w**2 * fen.dot(self.v1, self.v2)
+ ) * fen.dx)
+ normMat = fen.as_backend_type(normMatFen).mat()
+ return scsp.csr_matrix(normMat.getValuesCSR()[::-1],
+ shape = normMat.size)
+
+ def problemData(self):
+ """List of HF problem data."""
+ dataDict = {}
+ dataDict["forcingTerm"] = self.forcingTerm
+ dataDict["DirichletDatum"] = self.DirichletDatum
+ dataDict["NeumannDatum"] = self.NeumannDatum
+ dataDict["RobinDatum2"] = self.RobinDatum2
+ return dataDict
+
+ def setProblemData(self, dataDict:dict, k2:complex):
+ """
+ Set problem data.
+
+ Args:
+ dataDict: Dictionary of problem data.
+ k2: Parameter value.
+ """
+ self.wavenumber2 = k2
+ self.forcingTerm = dataDict["forcingTerm"]
+ self.DirichletDatum = dataDict["DirichletDatum"]
+ self.NeumannDatum = dataDict["NeumannDatum"]
+ self.RobinDatum2 = dataDict["RobinDatum2"]
+
+ def getLSBlocks(self) -> ("numpy 2D array" * 3, "numpy 1D array"):
+ """
+ Get blocks of linear system.
+
+ Returns:
+ Blocks of system (\sum_{j=0}^J k^j A_j = b)
+ """
+ self.assembleA()
+ self.assembleb()
+ return [self.AL + self.AR, - self.AM], [self.b]
+
+ def liftDirichletData(self):
+ """Lift Dirichlet datum."""
+ solLRe = fen.interpolate(self.u0Re, self.V)
+ solLIm = fen.interpolate(self.u0Im, self.V)
+ return np.array(solLRe.vector()) + 1.j * np.array(solLIm.vector())
+
+ def setupDerivativeComputation(self, j:int, up:"numpy 1D array",
+ upp : "numpy 1D array" = None):
+ """
+ Setup problem data to compute solution derivatives.
+
+ Args:
+ j: Derivative index.
+ up: Adjusted previous derivative.
+ upp: Adjusted pre-previous derivative.
+ """
+ upRe = fen.Function(self.V)
+ upIm = fen.Function(self.V)
+ upRe.vector()[:] = np.array(np.real(up), dtype = float)
+ upIm.vector()[:] = np.array(np.imag(up), dtype = float)
+ upRe, upIm = self.n2Re * upRe - self.n2Im * upIm,\
+ self.n2Re * upIm + self.n2Im * upRe
+ self.forcingTerm = [upRe, upIm]
+ self.DirichletDatum = fenZEROC
+ self.NeumannDatum = fenZEROC
+ self.RobinDatum2 = fenZEROC
+
+ @property
+ def wavenumber2(self):
+ """Value of k^2."""
+ return self._wavenumber2
+ @wavenumber2.setter
+ def wavenumber2(self, wavenumber2):
+ if hasattr(self, "A"): del self.A
+ if hasattr(self, "u"): del self.u
+ self._wavenumber2 = wavenumber2
+
+ @property
+ def diffusivity(self):
+ """Value of a. Its assignment changes aRe and aIm."""
+ return self._diffusivity
+ @diffusivity.setter
+ def diffusivity(self, diffusivity):
+ if hasattr(self, "A"): del self.A
+ if hasattr(self, "AL"): del self.AL
+ if hasattr(self, "u"): del self.u
+ self.aRe, self.aIm = CustomExpressionParser(diffusivity,
+ degree = self.FEDegree)
+ self._diffusivity = copy(diffusivity)
+
+ @property
+ def refractionIndex(self):
+ """Value of n. Its assignment changes nRe, nIm, n2Re and n2Im."""
+ return self._refractionIndex
+ @refractionIndex.setter
+ def refractionIndex(self, refractionIndex):
+ if hasattr(self, "A"): del self.A
+ if hasattr(self, "AM"): del self.AM
+ if hasattr(self, "u"): del self.u
+ self.nRe, self.nIm = CustomExpressionParser(refractionIndex,
+ degree = self.FEDegree)
+ self.n2Re = self.nRe*self.nRe - self.nIm*self.nIm
+ self.n2Im = 2 * self.nRe * self.nIm
+ self._refractionIndex = copy(refractionIndex)
+
+ @property
+ def forcingTerm(self):
+ """Value of f. Its assignment changes fRe and fIm."""
+ return self._forcingTerm
+ @forcingTerm.setter
+ def forcingTerm(self, forcingTerm):
+ if hasattr(self, "b"): del self.b
+ if hasattr(self, "bF"): del self.bF
+ if hasattr(self, "u"): del self.u
+ self.fRe, self.fIm = CustomExpressionParser(forcingTerm,
+ degree = self.FEDegree)
+ self._forcingTerm = copy(forcingTerm)
+
+ @property
+ def DirichletDatum(self):
+ """
+ Value of u0. Its assignment changes u0Re, u0Im, DirichletBCRe and
+ DirichletBCIm.
+ """
+ return self._DirichletDatum
+ @DirichletDatum.setter
+ def DirichletDatum(self, DirichletDatum):
+ if hasattr(self, "b"): del self.b
+ if hasattr(self, "u"): del self.u
+ self.u0Re, self.u0Im = CustomExpressionParser(DirichletDatum,
+ degree = self.FEDegree)
+ self.DirichletBCRe = fen.DirichletBC(self.V, self.u0Re,
+ self.DirichletBoundary)
+ self.DirichletBCIm = fen.DirichletBC(self.V, self.u0Im,
+ self.DirichletBoundary)
+ self.DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
+ self.DirichletBoundary)
+ self._DirichletDatum = copy(DirichletDatum)
+
+ @property
+ def NeumannDatum(self):
+ """Value of g1. Its assignment changes g1Re and g1Im."""
+ return self._NeumannDatum
+ @NeumannDatum.setter
+ def NeumannDatum(self, NeumannDatum):
+ if hasattr(self, "b"): del self.b
+ if hasattr(self, "bN"): del self.bN
+ if hasattr(self, "u"): del self.u
+ self.g1Re, self.g1Im = CustomExpressionParser(NeumannDatum,
+ degree = self.FEDegree)
+ self._NeumannDatum = copy(NeumannDatum)
+
+ @property
+ def RobinDatum(self):
+ """Value of h. Its assignment changes hRe and hIm."""
+ return self._RobinDatum
+ @RobinDatum.setter
+ def RobinDatum(self, RobinDatum):
+ if hasattr(self, "A"): del self.A
+ if hasattr(self, "AR"): del self.AR
+ if hasattr(self, "u"): del self.u
+ self.hRe, self.hIm = CustomExpressionParser(RobinDatum,
+ degree = self.FEDegree)
+ self._RobinDatum = copy(RobinDatum)
+
+ @property
+ def RobinDatum2(self):
+ """Value of g2. Its assignment changes g2Re and g2Im."""
+ return self._RobinDatum2
+ @RobinDatum2.setter
+ def RobinDatum2(self, RobinDatum2):
+ if hasattr(self, "b"): del self.b
+ if hasattr(self, "bR"): del self.bR
+ if hasattr(self, "u"): del self.u
+ self.g2Re, self.g2Im = CustomExpressionParser(RobinDatum2,
+ degree = self.FEDegree)
+ self._RobinDatum2 = copy(RobinDatum2)
+
+ def assembleA(self, Rfact : complex = 1.):
+ """Assemble matrix of linear system."""
+ if not hasattr(self, "A"):
+ if not hasattr(self, "AL"):
+ aLRe = self.aRe * fen.dot(fen.grad(self.v1),
+ fen.grad(self.v2)) * fen.dx
+ aLIm = self.aIm * fen.dot(fen.grad(self.v1),
+ fen.grad(self.v2)) * fen.dx
+ ALRe = fen.assemble(aLRe)
+ ALIm = fen.assemble(aLIm)
+ self.DirichletBC0.apply(ALRe)
+ self.DirichletBC0.zero(ALIm)
+ ALReMat = fen.as_backend_type(ALRe).mat()
+ ALRer, ALRec, ALRev = ALReMat.getValuesCSR()
+ ALImr, ALImc, ALImv = fen.as_backend_type(
+ ALIm).mat().getValuesCSR()
+ self.AL = (scsp.csr_matrix((ALRev, ALRec, ALRer),
+ shape = ALReMat.size)
+ + 1.j * scsp.csr_matrix((ALImv, ALImc, ALImr),
+ shape = ALReMat.size))
+ if not hasattr(self, "AM"):
+ aMRe = self.n2Re * fen.dot(self.v1, self.v2) * fen.dx
+ aMIm = self.n2Im * fen.dot(self.v1, self.v2) * fen.dx
+ AMRe = fen.assemble(aMRe)
+ AMIm = fen.assemble(aMIm)
+ self.DirichletBC0.zero(AMRe)
+ self.DirichletBC0.zero(AMIm)
+ AMRer, AMRec, AMRev = fen.as_backend_type(
+ AMRe).mat().getValuesCSR()
+ AMImr, AMImc, AMImv = fen.as_backend_type(
+ AMIm).mat().getValuesCSR()
+ self.AM = (scsp.csr_matrix((AMRev, AMRec, AMRer),
+ shape = self.AL.shape)
+ + 1.j * scsp.csr_matrix((AMImv, AMImc, AMImr),
+ shape = self.AL.shape))
+ if not hasattr(self, "AR"):
+ aRRe = self.hRe * fen.dot(self.v1, self.v2) * self.ds(1)
+ aRIm = self.hIm * fen.dot(self.v1, self.v2) * self.ds(1)
+ ARRe = fen.assemble(aRRe)
+ ARIm = fen.assemble(aRIm)
+ self.DirichletBC0.zero(ARRe)
+ self.DirichletBC0.zero(ARIm)
+ ARRer, ARRec, ARRev = fen.as_backend_type(
+ ARRe).mat().getValuesCSR()
+ ARImr, ARImc, ARImv = fen.as_backend_type(
+ ARIm).mat().getValuesCSR()
+ self.AR = (scsp.csr_matrix((ARRev, ARRec, ARRer),
+ shape = self.AL.shape)
+ + 1.j * scsp.csr_matrix((ARImv, ARImc, ARImr),
+ shape = self.AL.shape))
+
+ self.A = self.AL - self.wavenumber2 * self.AM + Rfact * self.AR
+ if hasattr(self, "invA"): del self.invA
+
+ def assembleb(self):
+ """Assemble RHS of linear system."""
+ if not hasattr(self, "b"):
+ if not hasattr(self, "bF"):
+ LFRe = fen.dot(self.fRe, self.v2) * fen.dx
+ LFIm = fen.dot(self.fIm, self.v2) * fen.dx
+ bFRe = fen.assemble(LFRe)
+ bFIm = fen.assemble(LFIm)
+ self.DirichletBCRe.apply(bFRe)
+ self.DirichletBCIm.apply(bFIm)
+ self.bF = np.array(bFRe.array()[:] + 1.j * bFIm.array()[:],
+ dtype = np.complex)
+ if not hasattr(self, "bN"):
+ LNRe = fen.dot(self.g1Re, self.v2) * self.ds(0)
+ LNIm = fen.dot(self.g1Im, self.v2) * self.ds(0)
+ bNRe = fen.assemble(LNRe)
+ bNIm = fen.assemble(LNIm)
+ self.DirichletBC0.apply(bNRe)
+ self.DirichletBC0.apply(bNIm)
+ self.bN = np.array(bNRe.array()[:] + 1.j * bNIm.array()[:],
+ dtype = np.complex)
+ if not hasattr(self, "bR"):
+ LRRe = fen.dot(self.g2Re, self.v2) * self.ds(1)
+ LRIm = fen.dot(self.g2Im, self.v2) * self.ds(1)
+ bRRe = fen.assemble(LRRe)
+ bRIm = fen.assemble(LRIm)
+ self.DirichletBC0.apply(bRRe)
+ self.DirichletBC0.apply(bRIm)
+ self.bR = np.array(bRRe.array()[:] + 1.j * bRIm.array()[:],
+ dtype = np.complex)
+
+ self.b = self.bF + self.bN + self.bR
+
+ def buildLU(self):
+ """Build LU decomposition of A."""
+ if not hasattr(self, "A"): self.assembleA()
+ if not hasattr(self, "invA"):
+ warnings.simplefilter("ignore", scsp.SparseEfficiencyWarning)
+ self.invA = spla.splu(self.A)
+ warnings.simplefilter("default", scsp.SparseEfficiencyWarning)
+
+ def solve(self, LU : bool = False)\
+ -> ("Fenics function", "Fenics function"):
+ """
+ Find solution of linear system.
+
+ Args:
+ LU(optional): Whether to use a LU solver for the system. Defaults
+ to False.
+
+ Returns:
+ Real and imaginary parts of solution.
+ """
+ if not hasattr(self, "u"):
+ if not hasattr(self, "A"): self.assembleA()
+ if not hasattr(self, "b"): self.assembleb()
+ if LU:
+ self.buildLU()
+ self.u = self.invA.solve(self.b)
+ else:
+ self.u = spla.spsolve(self.A, self.b)
+ self.__solved = True
+ return self.u
+
+
+class FenicsHelmholtzScatteringEngine(FenicsHelmholtzEngine):
+ """
+ Fenics-based solver for Helmholtz scattering problems with Robin ABC.
+ - \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega
+ u = u0 on \Gamma_D
+ \partial_nu = g1 on \Gamma_N
+ \partial_nu - i k u = g2 on \Gamma_R
+
+ Args:
+ mesh: Domain of Helmholtz problem.
+ FEDegree: FE degree.
+ wavenumber: Value of k.
+ diffusivity(optional): Value of a. Defaults to identically 1.
+ refractionIndex(optional): Value of n. Defaults to identically 1.
+ forcingTerm(optional): Value of f. Defaults to identically 0.
+ DirichletBoundary(optional): Function flagging Dirichlet boundary as
+ True, in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are
+ accepted. Defaults to False everywhere.
+ NeumannBoundary(optional): Function flagging Neumann boundary as True,
+ in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
+ Defaults to False everywhere.
+ RobinBoundary(optional): Function flagging Robin boundary as True, in
+ Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
+ Defaults to False everywhere.
+ DirichletDatum(optional): Value of u0. Defaults to identically 0.
+ NeumannDatum(optional): Value of g1. Defaults to identically 0.
+ RobinDatum2(optional): Value of g2. Defaults to identically 0.
+
+ Attributes:
+ mesh: Domain of Helmholtz problem.
+ FEDegree: FE degree.
+ wavenumber: Copy of processed wavenumber parameter.
+ diffusivity: Copy of processed diffusivity parameter.
+ refractionIndex: Copy of processed refractionIndex parameter.
+ forcingTerm: Copy of processed forcingTerm parameter.
+ DirichletBoundary: Copy of processed DirichletBoundary parameter.
+ NeumannBoundary: Copy of processed NeumannBoundary parameter.
+ RobinBoundary: Copy of processed RobinBoundary parameter.
+ DirichletDatum: Copy of processed DirichletDatum parameter.
+ NeumannDatum: Copy of processed NeumannDatum parameter.
+ RobinDatum: Value of h, i.e. (- i * k).
+ RobinDatum2: Copy of processed RobinDatum2 parameter.
+ V: Real FE space.
+ u: Helmholtz problem solution as numpy complex vector.
+ v1: Generic trial functions for variational form evaluation.
+ v2: Generic test functions for variational form evaluation.
+ ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
+ nRe,nIm: Real and imaginary parts of n.
+ n2Re,n2Im: Real and imaginary parts of n^2.
+ fRe,fIm: Real and imaginary parts of f.
+ u0Re,u0Im: Real and imaginary parts of u0.
+ g1Re,g1Im: Real and imaginary parts of g1.
+ hRe,hIm: Real and imaginary parts of h.
+ g2Re,g2Im: Real and imaginary parts of g2.
+ DirichletBCRe,DirichletBCIm: Fenics BC manager for real and imaginary
+ parts of Dirichlet data.
+ A*: Scipy sparse array representation (in CSC format) of A*.
+ b*Re,b*Im: Numpy array representation of real and imaginary parts of
+ b*.
+ solRe,solIm: Real and imaginary parts of Helmholtz problem solution.
+ invA: Numpy LU solver for A.
+ """
+
+ def __init__(self, mesh:"Fenics mesh",
+ FEDegree:int,
+ wavenumber:"custom complex",
+ diffusivity : "custom expression" = 1,
+ refractionIndex : "custom expression" = 1,
+ forcingTerm : "custom expression" = fenZEROC,
+ DirichletBoundary : "bool function or str" = None,
+ NeumannBoundary : "bool function or str" = None,
+ RobinBoundary : "bool function or str" = None,
+ DirichletDatum : "custom expression" = fenZEROC,
+ NeumannDatum : "custom expression" = fenZEROC,
+ RobinDatum2 : "custom expression" = fenZEROC):
+ self.wavenumber = wavenumber
+ FenicsHelmholtzEngine.__init__(self, mesh = mesh,
+ FEDegree = FEDegree,
+ wavenumber = wavenumber,
+ diffusivity = diffusivity,
+ refractionIndex = refractionIndex,
+ forcingTerm = forcingTerm,
+ DirichletBoundary = DirichletBoundary,
+ NeumannBoundary = NeumannBoundary,
+ RobinBoundary = RobinBoundary,
+ DirichletDatum = DirichletDatum,
+ NeumannDatum = NeumannDatum,
+ RobinDatum = -1.j * wavenumber,
+ RobinDatum2 = RobinDatum2)
+
+ def setProblemData(self, dataDict:dict, k:complex):
+ """
+ Set problem data.
+
+ Args:
+ dataDict: Dictionary of problem data.
+ k: Parameter value.
+ """
+ self.wavenumber = k
+ self.forcingTerm = dataDict["forcingTerm"]
+ self.DirichletDatum = dataDict["DirichletDatum"]
+ self.NeumannDatum = dataDict["NeumannDatum"]
+ self.RobinDatum2 = dataDict["RobinDatum2"]
+
+ def getLSBlocks(self) -> ("numpy 2D array" * 3, "numpy 1D array"):
+ """
+ Get blocks of linear system.
+
+ Returns:
+ Blocks of system (\sum_{j=0}^J k^j A_j = b)
+ """
+ self.assembleA()
+ self.assembleb()
+ return [self.AL, - 1.j * self.AR, - self.AM], [self.b]
+
+ def setupDerivativeComputation(self, j:int, up:"numpy 1D array",
+ upp : "numpy 1D array"):
+ """
+ Setup problem data to compute solution derivatives.
+
+ Args:
+ j: Derivative index.
+ up: Adjusted previous derivative.
+ upp: Adjusted pre-previous derivative.
+ """
+ upRe = fen.Function(self.V)
+ upIm = fen.Function(self.V)
+ ftRe = fen.Function(self.V)
+ ftIm = fen.Function(self.V)
+ upRe.vector()[:] = np.array(np.real(up), dtype=float)
+ upIm.vector()[:] = np.array(np.imag(up), dtype=float)
+ ft0 = 2 * self.wavenumber * up + upp
+ ftRe.vector()[:] = np.array(np.real(ft0), dtype=float)
+ ftIm.vector()[:] = np.array(np.imag(ft0), dtype=float)
+ ftRe, ftIm = self.n2Re * ftRe - self.n2Im * ftIm,\
+ self.n2Re * ftIm + self.n2Im * ftRe
+ self.forcingTerm = [ftRe, ftIm]
+ self.DirichletDatum = fenZEROC
+ self.NeumannDatum = fenZEROC
+ self.RobinDatum2 = [-upIm, upRe]
+
+ @property
+ def wavenumber2(self):
+ """Value of k^2."""
+ return self._wavenumber2
+ @wavenumber2.setter
+ def wavenumber2(self, wavenumber2):
+ FenicsHelmholtzEngine.wavenumber2.fset(self, wavenumber2)
+ self._wavenumber = self.wavenumber2 ** .5
+ if (not hasattr(self, "h")
+ or not np.isclose(self.h, -1.j * self.wavenumber, 1e-14)):
+ self.RobinDatum = -1.j * self.wavenumber
+
+ @property
+ def wavenumber(self):
+ """Value of k."""
+ return self._wavenumber
+ @wavenumber.setter
+ def wavenumber(self, wavenumber):
+ self.wavenumber2 = wavenumber ** 2.
+
+ @property
+ def RobinDatum(self):
+ """
+ Value of h. Its assignment changes hRe and hIm (and maybe kRe, kIm, k
+ and z).
+ """
+ return super().RobinDatum
+ @RobinDatum.setter
+ def RobinDatum(self, RobinDatum):
+ if hasattr(self, "A"): del self.A
+ if hasattr(self, "solRe"): del self.solRe, self.solIm
+ self.h = RobinDatum
+ self._RobinDatum = copy(RobinDatum)
+ if (not hasattr(self, "wavenumber")
+ or not np.isclose(self.h, -1.j * self.wavenumber, 1e-14)):
+ self.wavenumber = 1.j * self.h
+
+ def assembleA(self):
+ """Assemble matrix of linear system."""
+ if not hasattr(self, "A"):
+ if not hasattr(self, "AR"):
+ aR = fen.dot(self.v1, self.v2) * self.ds(1)
+ AR = fen.assemble(aR)
+ self.DirichletBC0.zero(AR)
+ AR_mat = fen.as_backend_type(AR).mat().transpose()
+ self.AR = scsp.csc_matrix(AR_mat.getValuesCSR()[::-1],
+ shape = AR_mat.size)
+ FenicsHelmholtzEngine.assembleA(self, self.h)
+
+
+class FenicsHelmholtzScatteringAugmentedEngine(FenicsHelmholtzScatteringEngine):
+ """
+ Fenics-based solver for Helmholtz scattering problems with Robin ABC.
+ - \nabla \cdot (a \nabla u) - k^2 * n**2 * u = f in \Omega
+ u = u0 on \Gamma_D
+ \partial_nu = g1 on \Gamma_N
+ \partial_nu - i k u = g2 on \Gamma_R
+ Linear dependence on k is achieved by introducing an additional variable.
+
+ Args:
+ mesh: Domain of Helmholtz problem.
+ FEDegree: FE degree.
+ wavenumber: Value of k.
+ diffusivity(optional): Value of a. Defaults to identically 1.
+ refractionIndex(optional): Value of n. Defaults to identically 1.
+ forcingTerm(optional): Value of f. Defaults to identically 0.
+ DirichletBoundary(optional): Function flagging Dirichlet boundary as
+ True, in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are
+ accepted. Defaults to False everywhere.
+ NeumannBoundary(optional): Function flagging Neumann boundary as True,
+ in Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
+ Defaults to False everywhere.
+ RobinBoundary(optional): Function flagging Robin boundary as True, in
+ Fenics format. Keywords 'ALL', 'NONE' and 'REST' are accepted.
+ Defaults to False everywhere.
+ DirichletDatum(optional): Value of u0. Defaults to identically 0.
+ NeumannDatum(optional): Value of g1. Defaults to identically 0.
+ RobinDatum2(optional): Value of g2. Defaults to identically 0.
+ constraintType(optional): Type of augmentation. Keywords 'IDENTITY' and
+ 'MASS' are accepted. Defaults to 'IDENTITY'.
+
+ Attributes:
+ mesh: Domain of Helmholtz problem.
+ FEDegree: FE degree.
+ wavenumber: Copy of processed wavenumber parameter.
+ diffusivity: Copy of processed diffusivity parameter.
+ refractionIndex: Copy of processed refractionIndex parameter.
+ forcingTerm: Copy of processed forcingTerm parameter.
+ DirichletBoundary: Copy of processed DirichletBoundary parameter.
+ NeumannBoundary: Copy of processed NeumannBoundary parameter.
+ RobinBoundary: Copy of processed RobinBoundary parameter.
+ DirichletDatum: Copy of processed DirichletDatum parameter.
+ NeumannDatum: Copy of processed NeumannDatum parameter.
+ RobinDatum: Value of h, i.e. (- i * k).
+ RobinDatum2: Copy of processed RobinDatum2 parameter.
+ constraintType: Type of augmentation.
+ V: Real FE space.
+ u: Helmholtz problem solution as numpy complex vector.
+ v1: Generic trial functions for variational form evaluation.
+ v2: Generic test functions for variational form evaluation.
+ ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
+ nRe,nIm: Real and imaginary parts of n.
+ n2Re,n2Im: Real and imaginary parts of n^2.
+ fRe,fIm: Real and imaginary parts of f.
+ u0Re,u0Im: Real and imaginary parts of u0.
+ g1Re,g1Im: Real and imaginary parts of g1.
+ hRe,hIm: Real and imaginary parts of h.
+ g2Re,g2Im: Real and imaginary parts of g2.
+ DirichletBCRe,DirichletBCIm: Fenics BC manager for real and imaginary
+ parts of Dirichlet data.
+ A*: Scipy sparse array representation (in CSC format) of A*.
+ b*Re,b*Im: Numpy array representation of real and imaginary parts of
+ b*.
+ solRe,solIm: Real and imaginary parts of Helmholtz problem solution.
+ invA: Numpy LU solver for A.
+ """
+
+ def __init__(self, mesh:"Fenics mesh",
+ FEDegree:int,
+ wavenumber:"custom complex",
+ diffusivity : "custom expression" = 1,
+ refractionIndex : "custom expression" = 1,
+ forcingTerm : "custom expression" = fenZEROC,
+ DirichletBoundary : "bool function or str" = None,
+ NeumannBoundary : "bool function or str" = None,
+ RobinBoundary : "bool function or str" = None,
+ DirichletDatum : "custom expression" = fenZEROC,
+ NeumannDatum : "custom expression" = fenZEROC,
+ RobinDatum2 : "custom expression" = fenZEROC,
+ constraintType : str = "IDENTITY"):
+ self.wavenumber = wavenumber
+ self.constraintType = constraintType
+ FenicsHelmholtzScatteringEngine.__init__(self, mesh = mesh,
+ FEDegree = FEDegree,
+ wavenumber = wavenumber,
+ diffusivity = diffusivity,
+ refractionIndex = refractionIndex,
+ forcingTerm = forcingTerm,
+ DirichletBoundary = DirichletBoundary,
+ NeumannBoundary = NeumannBoundary,
+ RobinBoundary = RobinBoundary,
+ DirichletDatum = DirichletDatum,
+ NeumannDatum = NeumannDatum,
+ RobinDatum2 = RobinDatum2)
+
+ def energyNormMatrix(self, w:float) -> "CSC sparse matrix":
+ """
+ Sparse matrix (in CSR format) assoociated to w-weighted H10 norm.
+
+ Args:
+ w: Weight.
+
+ Returns:
+ Sparse matrix in CSR format.
+ """
+ normMatFen = FenicsHelmholtzEngine.energyNormMatrix(self, w)
+ return scsp.block_diag((normMatFen, 1/w * normMatFen))
+
+ def liftDirichletData(self):
+ """Lift Dirichlet datum."""
+ lift1 = FenicsHelmholtzScatteringEngine.liftDirichletData(self)
+ return np.pad(lift1, (0, len(lift1)), 'constant')
+
+ @property
+ def constraintType(self):
+ """Value of constraintType."""
+ return self._constraintType
+ @constraintType.setter
+ def constraintType(self, constraintType):
+ try:
+ constraintType = constraintType.upper().strip().replace(" ","")
+ if constraintType not in ["IDENTITY", "MASS"]: raise
+ self._constraintType = constraintType
+ except:
+ warnings.warn(("Prescribed constraintType not recognized. "
+ "Overriding to 'KRYLOV'."))
+ self.constraintType = "IDENTITY"
+
+ def getLSBlocks(self) -> (list, list):
+ """
+ Get blocks of linear system.
+
+ Returns:
+ Blocks of system (\sum_{j=0}^J k^j A_j = b)
+ """
+ self.assembleA()
+ self.assembleb()
+ return [self.A0, - self.A1], [self.b]
+
+ def setupDerivativeComputation(self, j:int, up:"numpy 1D array",
+ upp : "numpy 1D array"):
+ """
+ Setup problem data to compute solution derivatives.
+
+ Args:
+ j: Derivative index.
+ up: Adjusted previous derivative.
+ upp: Adjusted pre-previous derivative.
+ """
+ up1Re = fen.Function(self.V)
+ up1Im = fen.Function(self.V)
+ up2Re = fen.Function(self.V)
+ up2Im = fen.Function(self.V)
+ lup = int(len(up) / 2)
+ up1Re.vector()[:] = np.array(np.real(up[: lup]), dtype=float)
+ up1Im.vector()[:] = np.array(np.imag(up[: lup]), dtype=float)
+ up2Re.vector()[:] = np.array(np.real(up[lup :]), dtype=float)
+ up2Im.vector()[:] = np.array(np.imag(up[lup :]), dtype=float)
+ if self.constraintType == "IDENTITY":
+ b2 = up[: lup]
+ else:
+ self.forcingTerm = [up1Re, up1Im]
+ self.DirichletDatum = fenZEROC
+ self.NeumannDatum = fenZEROC
+ self.RobinDatum2 = fenZEROC
+ self.assembleb()
+ b2 = self.b[lup :]
+ self.forcingTerm = [up2Re, up2Im]
+ self.DirichletDatum = fenZEROC
+ self.NeumannDatum = fenZEROC
+ self.RobinDatum2 = [-up1Im, up1Re]
+ self.assembleb()
+ self.b[lup :] = b2
+
+ def assembleA(self):
+ """Assemble matrix of linear system."""
+ if not hasattr(self, "A"):
+ if hasattr(self, "A0") and not hasattr(self, "AL"):
+ del self.A0
+ if (hasattr(self, "A1")
+ and not (hasattr(self, "AM") and hasattr(self, "AR"))):
+ del self.A1
+ FenicsHelmholtzScatteringEngine.assembleA(self)
+ if self.constraintType == "IDENTITY":
+ block1 = scsp.eye(self.AL.shape[0])
+ block2 = block1
+ else: #MASS
+ block1 = copy(self.AM)
+ I = list(self.DirichletBC0.get_boundary_values().keys())
+ warnings.simplefilter("ignore", scsp.SparseEfficiencyWarning)
+ block1[I, I] = 1.
+ warnings.simplefilter("default", scsp.SparseEfficiencyWarning)
+ block2 = self.AM
+ self.A0 = scsp.block_diag((self.AL, block1))
+ if not hasattr(self, "A1"):
+ self.A1 = scsp.bmat([[1.j * self.AR, self.AM], [block2, None]])
+
+ self.A = self.A0 - self.wavenumber * self.A1
+ if hasattr(self, "invA"): del self.invA
+
+ def assembleb(self):
+ """Assemble RHS of linear system."""
+ if not hasattr(self, "b"):
+ FenicsHelmholtzScatteringEngine.assembleb(self)
+ self.b = np.pad(self.b, (0, self.b.shape[0]), 'constant')
diff --git a/main/ROMApproximant.py b/main/ROMApproximant.py
new file mode 100644
index 0000000..c0069d3
--- /dev/null
+++ b/main/ROMApproximant.py
@@ -0,0 +1,307 @@
+#!/usr/bin/python
+
+from abc import abstractmethod
+import numpy as np
+import utilities
+
+class ROMApproximant:
+ """
+ ROM approximant computation for parametric problems.
+
+ Args:
+ HFEngine: HF problem solver. Should have members:
+ - energyNormMatrix: assemble sparse matrix (in CSC format)
+ associated to weighted H10 norm;
+ - problemData: list of HF problem data (excluding k);
+ - setProblemData: set HF problem data and k to prescribed values;
+ - solve: get HF solution as complex numpy vector.
+ HSEngine: Hilbert space general purpose engine. Should have members:
+ - norm: compute Hilbert norm of expression represented as complex
+ numpy vector;
+ - plot: plot Hilbert expression represented as complex numpy vector.
+ k0(optional): Default parameter. Defaults to 0.
+ w(optional): Weight for norm computation. If None, set to Re(k0).
+ Defaults to None.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True.
+ Defaults to empty dict.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ HSEngine: Hilbert space general purpose engine.
+ k0: Default parameter.
+ w: Weight for norm computation.
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots.
+ initialHFData: HF problem initial data.
+ energyNormMatrix: Sparse matrix (in CSC format) assoociated to
+ w-weighted H10 norm.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ """
+
+ def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
+ k0 : complex = 0, w : float = None,
+ approxParameters : dict = {}):
+ self.HFEngine = HFEngine
+ self.HSEngine = HSEngine
+ self.initialHFData = self.HFEngine.problemData()
+
+ self.k0 = k0
+ if w is None:
+ w = np.real(self.k0)
+ self.w = w
+ self.approxParameters = approxParameters
+ self.energyNormMatrix = self.HFEngine.energyNormMatrix(self.w)
+
+ def name(self) -> str:
+ """Approximant label."""
+ return self.__class__.__name__
+
+ def parameterList(self) -> list:
+ """
+ List containing keys which are allowed in approxParameters.
+
+ Returns:
+ List of strings.
+ """
+ return ["POD"]
+
+ @property
+ def approxParameters(self):
+ """Value of approximant parameters."""
+ return self._approxParameters
+ @approxParameters.setter
+ def approxParameters(self, approxParams):
+ if not hasattr(self, "approxParameters"):
+ self._approxParameters = {}
+ approxParameters = utilities.purgeDict(approxParams,
+ ROMApproximant.parameterList(self),
+ dictname = "approxParameters")
+ keyList = list(approxParameters.keys())
+ if "POD" in keyList:
+ self.POD = approxParameters["POD"]
+ elif hasattr(self, "POD"):
+ self.POD = self.POD
+ else:
+ self.POD = True
+
+ @property
+ def POD(self):
+ """Value of POD."""
+ return self._POD
+ @POD.setter
+ def POD(self, POD:bool):
+ if hasattr(self, "POD"): PODold = self.POD
+ else: PODold = -1
+ self._POD = POD
+ self._approxParameters["POD"] = self.POD
+ if PODold != self.POD:
+ self.resetSamples()
+
+ def solveHF(self, k : complex = None):
+ """
+ Find high fidelity solution with original parameters and arbitrary
+ wavenumber.
+
+ Args:
+ k: Target wavenumber.
+ """
+ if k is None: k = self.k0
+ if (not hasattr(self, "lastSolvedHF")
+ or not np.isclose(self.lastSolvedHF, k)):
+ self.HFEngine.setProblemData(self.initialHFData, k)
+ self.uHF = self.HFEngine.solve(False)
+
+ @abstractmethod
+ def resetSamples(self):
+ """Reset samples. (ABSTRACT)"""
+ pass
+
+ @abstractmethod
+ def setupApprox(self):
+ """
+ Setup approximant. (ABSTRACT)
+
+ Any specialization should include something like
+ self.computeDerivatives()
+ if not self.checkComputedApprox():
+ ...
+ self.lastApproxParameters = copy(self.approxParameters)
+ """
+ pass
+
+ def checkComputedApprox(self) -> bool:
+ """
+ Check if setup of new approximant is not needed.
+
+ Returns:
+ True if new setup is not needed. False otherwise.
+ """
+ return (hasattr(self, "lastApproxParameters")
+ and self.approxParameters == self.lastApproxParameters)
+
+ @abstractmethod
+ def evalApprox(self, k:complex) -> ("numpy 1D array"):
+ """
+ Evaluate approximant at arbitrary wavenumber. (ABSTRACT)
+
+ Any specialization should include something like
+ self.setupApprox()
+ self.uApp = ...
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ Approximant as numpy complex vector.
+ """
+ pass
+
+ @abstractmethod
+ def getPoles(self, centered : bool = False) -> "numpy 1D array":
+ """
+ Obtain approximant poles.
+
+ Args:
+ centered(optional): Whether to return pole positions relative to
+ approximation center. Defaults to False.
+
+ Returns:
+ Numpy complex vector of poles.
+ """
+ pass
+
+ def HFNorm(self, k:complex, normType : "number or str" = None) -> float:
+ """
+ Compute norm of HF solution at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+ normType(optional): Target norm identifier. If number, target norm
+ is weighted H^1 norm with given weight. If string, must be
+ recognizable by Fenics norm command. If None, set to w.
+ Defaults to None.
+
+ Returns:
+ Target norm of HFsolution.
+ """
+ self.solveHF(k)
+ if normType is None: normType = self.w
+ return self.HSEngine.norm(self.uHF, normType)
+
+ def approxNorm(self, k:complex, normType : "number or str" = None)-> float:
+ """
+ Compute norm of (translated) approximant at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+ normType(optional): Target norm identifier. If number, target norm
+ is weighted H^1 norm with given weight. If string, must be
+ recognizable by Fenics norm command. If None, set to w.
+ Defaults to None.
+
+ Returns:
+ Target norm of approximant.
+ """
+ self.evalApprox(k)
+ if normType is None: normType = self.w
+ return self.HSEngine.norm(self.uApp, normType)
+
+ def approxError(self, k:complex, normType : "number or str" = None)\
+ -> float:
+ """
+ Compute norm of approximant at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+ normType(optional): Target norm identifier. If number, target norm
+ is weighted H^1 norm with given weight. If string, must be
+ recognizable by Fenics norm command. If None, set to w.
+ Defaults to None.
+
+ Returns:
+ Target norm of (approximant - HFsolution).
+ """
+ self.evalApprox(k)
+ self.solveHF(k)
+ self.evalApprox(k)
+ if normType is None: normType = self.w
+ return self.HSEngine.norm(self.uApp - self.uHF, normType)
+
+ def getHF(self, k:complex) -> "numpy 1D array":
+ """
+ Get HF solution at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ HFsolution as numpy complex vector.
+ """
+ self.solveHF(k)
+ return self.uHF
+
+ def getApp(self, k:complex) -> "numpy 1D array":
+ """
+ Get approximant at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ Approximant as numpy complex vector.
+ """
+ self.evalApprox(k)
+ return self.uApp
+
+ def plotHF(self, k:complex, name : str = "u", **figspecs):
+ """
+ Do some nice plots of the HF solution at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+ name(optional): Name to be shown as title of the plots. Defaults to
+ 'u'.
+ figspecs(optional key args): Optional arguments for matplotlib
+ figure creation.
+ """
+ self.solveHF(k)
+ self.HSEngine.plot(self.uHF, name, **figspecs)
+
+ def plotApp(self, k:complex, name : str = "u", **figspecs):
+ """
+ Do some nice plots of approximant at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+ name(optional): Name to be shown as title of the plots. Defaults to
+ 'u'.
+ figspecs(optional key args): Optional arguments for matplotlib
+ figure creation.
+ """
+ self.evalApprox(k)
+ self.HSEngine.plot(self.uApp, name, **figspecs)
+
+ def plotErr(self, k:complex, name : str = "u", **figspecs):
+ """
+ Do some nice plots of approximation error at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+ name(optional): Name to be shown as title of the plots. Defaults to
+ 'u'.
+ figspecs(optional key args): Optional arguments for matplotlib
+ figure creation.
+ """
+ self.evalApprox(k)
+ self.solveHF(k)
+ self.HSEngine.plot(self.uApp - self.uHF, name, **figspecs)
+
diff --git a/main/ROMApproximantLagrange.py b/main/ROMApproximantLagrange.py
new file mode 100644
index 0000000..b1da41f
--- /dev/null
+++ b/main/ROMApproximantLagrange.py
@@ -0,0 +1,193 @@
+#!/usr/bin/python
+
+import numpy as np
+import utilities
+from ROMApproximant import ROMApproximant
+
+
+class ROMApproximantLagrange(ROMApproximant):
+ """
+ ROM Lagrange interpolant computation for parametric problems.
+
+ Args:
+ HFEngine: HF problem solver. Should have members:
+ - energyNormMatrix: assemble sparse matrix (in CSC format)
+ associated to weighted H10 norm;
+ - problemData: list of HF problem data (excluding k);
+ - setProblemData: set HF problem data and k to prescribed values;
+ - solve: get HF solution as complex numpy vector.
+ HSEngine: Hilbert space general purpose engine. Should have members:
+ - norm: compute Hilbert norm of expression represented as complex
+ numpy vector;
+ - plot: plot Hilbert expression represented as complex numpy vector.
+ ks(optional): Array of snapshot parameters. Defaults to np.array([0]).
+ w(optional): Weight for norm computation. If None, set to
+ Re(np.mean(ks)). Defaults to None.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True;
+ - 'S': total number of samples current approximant relies upon.
+ Defaults to empty dict.
+ plotSnapshots(optional): Whether to plot snapshots of the Helmholtz
+ solution map. Defaults to False.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ HSEngine: Hilbert space general purpose engine.
+ solSnapshots: Array whose columns are FE dofs of snapshots of solution
+ map.
+ k0: Default parameter.
+ ks: Array of snapshot parameters.
+ w: Weight for norm computation.
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots;
+ - 'S': total number of snapshots current approximant relies upon.
+ S: Number of solution snapshots over which current approximant is
+ based upon.
+ plotSnapshots: Whether to plot snapshots of the Helmholtz solution map.
+ POD: Whether to compute POD of snapshots.
+ energyNormMatrix: Sparse matrix (in CSC format) assoociated to
+ w-weighted H10 norm.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ """
+
+ extraApproxParameters = ["S"]
+
+ def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
+ ks : "Numpy 1D array" = np.array([0]), w : float = None,
+ approxParameters : dict = {}, plotSnapshots : bool = False):
+ self.ks = ks
+ ROMApproximant.__init__(self, HFEngine, HSEngine,
+ np.mean(ks), w, approxParameters)
+ if w is None:
+ w = np.mean(self.ksRe)
+ self.w = w
+ self.plotSnapshots = plotSnapshots
+
+ def parameterList(self) -> list:
+ """
+ List containing keys which are allowed in approxParameters.
+
+ Returns:
+ List of strings.
+ """
+ return (ROMApproximant.parameterList(self)
+ + ROMApproximantLagrange.extraApproxParameters)
+
+ @property
+ def ks(self):
+ """Value of ks. Its assignment may reset snapshots."""
+ return self._ks
+ @ks.setter
+ def ks(self, ks):
+ if hasattr(self, 'ks'):
+ ksOld = self.ks
+ else:
+ ksOld = None
+
+ self._ks = np.array(ks)
+ if (ksOld is None or len(self.ks) != len(ksOld)
+ or not np.allclose(self.ks, ksOld, 1e-14)):
+ self.resetSamples()
+
+ @property
+ def ksRe(self):
+ """Real part of ks."""
+ return np.real(self.ks)
+
+ @property
+ def ksIm(self):
+ """Imaginary part of ks."""
+ return np.imag(self.ks)
+
+ @property
+ def zs(self):
+ """Square of ks."""
+ return np.power(self.ks, 2.)
+
+ @property
+ def approxParameters(self):
+ """Value of approximant parameters. Its assignment may change S."""
+ return self._approxParameters
+ @approxParameters.setter
+ def approxParameters(self, approxParams):
+ if not hasattr(self, "approxParameters"):
+ self._approxParameters = {}
+ approxParameters = utilities.purgeDict(approxParams,
+ ROMApproximantLagrange.parameterList(self),
+ dictname = "approxParameters")
+ keyList = list(approxParameters.keys())
+ approxParametersCopy = utilities.purgeDict(approxParameters,
+ ROMApproximantLagrange.extraApproxParameters,
+ True, True)
+ ROMApproximant.approxParameters.fset(self, approxParametersCopy)
+ if "S" in keyList:
+ self.S = approxParameters["S"]
+ elif hasattr(self, "S"):
+ self.S = self.S
+ else:
+ self.S = 2
+
+ @property
+ def S(self):
+ """Value of S."""
+ return self._S
+ @S.setter
+ def S(self, S):
+ if S <= 0: raise ArithmeticError("S must be positive.")
+ if hasattr(self, "S"): Sold = self.S
+ else: Sold = -1
+ self._S = S
+ self._approxParameters["S"] = self.S
+ if Sold != self.S:
+ self.resetSamples()
+
+ def resetSamples(self):
+ """Reset samples. (ABSTRACT)"""
+ self.solSnapshots = None
+ self.RPOD = None
+
+ def computeSnapshots(self):
+ """
+ Compute snapshots of solution map.
+ """
+ if self.solSnapshots is None:
+ for j, k in enumerate(self.ks):
+ self.solveHF(k)
+ if self.plotSnapshots:
+ self.HSEngine.plot(self.uHF, name = "u({:.4f})".format(k))
+ self.manageSnapshots(self.uHF, j)
+
+ def manageSnapshots(self, u:"numpy 1D array", pos:int):
+ """
+ Store snapshots of solution map.
+
+ Args:
+ u: solution derivative as numpy complex vector;
+ pos: Derivative index.
+ """
+ if pos == 0:
+ self.solSnapshots = np.empty((u.shape[0], self.S),
+ dtype = np.complex)
+ self.solSnapshots[:, pos] = u
+ if self.POD:
+ if pos == 0:
+ self.RPOD = np.eye(self.S, dtype = np.complex)
+ beta = 1
+ for j in range(2):
+ nu = self.solSnapshots[:, pos].conj().dot(
+ self.energyNormMatrix.dot(self.solSnapshots[:, : pos])).conj()
+ self.RPOD[: pos, pos] = self.RPOD[: pos, pos] + beta * nu
+ eta = (self.solSnapshots[:, pos]
+ - self.solSnapshots[:, : pos].dot(nu))
+ beta = eta.conj().dot(self.energyNormMatrix.dot(eta))**.5
+ self.solSnapshots[:, pos] = eta / beta
+ self.RPOD[pos, pos] = beta * self.RPOD[pos, pos]
+
diff --git a/main/ROMApproximantLagrangePade.py b/main/ROMApproximantLagrangePade.py
new file mode 100644
index 0000000..457c2c1
--- /dev/null
+++ b/main/ROMApproximantLagrangePade.py
@@ -0,0 +1,453 @@
+#!/usr/bin/python
+
+from __future__ import print_function
+from copy import copy
+import warnings
+import numpy as np
+import utilities
+from ROMApproximantLagrange import ROMApproximantLagrange
+
+PI = np.pi
+
+class ROMApproximantLagrangePade(ROMApproximantLagrange):
+ """
+ ROM Lagrange Pade' interpolant computation for parametric problems.
+
+ Args:
+ HFEngine: HF problem solver. Should have members:
+ - energyNormMatrix: assemble sparse matrix (in CSC format)
+ associated to weighted H10 norm;
+ - problemData: list of HF problem data (excluding k);
+ - setProblemData: set HF problem data and k to prescribed values;
+ - solve: get HF solution as complex numpy vector.
+ HSEngine: Hilbert space general purpose engine. Should have members:
+ - norm: compute Hilbert norm of expression represented as complex
+ numpy vector;
+ - plot: plot Hilbert expression represented as complex numpy vector.
+ ks(optional): Array of snapshot parameters. Defaults to np.array([0]).
+ ws(optional): Array of snapshots weights. Defaults to uniform = 1.
+ w(optional): Weight for norm computation. If None, set to
+ Re(np.mean(ks)). Defaults to None.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True;
+ - 'S': total number of samples current approximant relies upon;
+ defaults to 2;
+ - 'M': degree of Pade' interpolant numerator; defaults to 0;
+ - 'N': degree of Pade' interpolant denominator; defaults to 0.
+ - 'polyBasis': label of polynomial basis for LS problem; available
+ values are:
+ - 'CHEBYSHEV': Chebyshev nodes and weights;
+ - 'GAUSSLEGENDRE': Gauss-Legendre nodes and weights;
+ defaults to 'CHEBYSHEV'.
+ Defaults to empty dict.
+ plotSnapshots(optional): Whether to plot snapshots of the Helmholtz
+ solution map. Defaults to False.
+ verboseRobust(optional): Verbosity level for robustness-related
+ parameter modifications. Defaults to 1.
+ cleanupParameters(optional): Parameter values for pole cleanup.
+ Available fields are:
+ - 'boolCondition'(bool function handle): if evaluation on pole
+ returns False, pole is removed; defaults to always True;
+ - 'residueCheck'(bool): whether to apply residue check; defaults to
+ False;
+ - 'residueNPoints'(int): number of sample points to be used in
+ residue estimation; defaults to 2;
+ - 'residueRadius'(float): sample radius to be used in residue
+ estimation by Cauchy formula; defaults to 1e-5;
+ - 'residueTol'(float): target tolerance to be used in residue
+ estimation; defaults to 1e-4.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ HSEngine: Hilbert space general purpose engine.
+ solSnapshots: Array whose columns are FE dofs of snapshots of solution
+ map.
+ k0: Default parameter.
+ ks: Array of snapshot parameters.
+ ws: Array of snapshots weights.
+ w: Weight for norm computation.
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots;
+ - 'S': total number of samples current approximant relies upon;
+ - 'M': degree of Pade' interpolant numerator;
+ - 'N': degree of Pade' interpolant denominator;
+ - 'polyBasis': label of polynomial basis for LS problem.
+ M: Numerator degree of approximant.
+ N: Denominator degree of approximant.
+ S: Number of solution snapshots over which current approximant is
+ based upon.
+ polyBasis: Polynomial basis for LS problem.
+ plotSnapshots: Whether to plot snapshots of the Helmholtz solution map.
+ POD: Whether to compute POD of snapshots.
+ verboseRobust: Verbosity level for robustness-related parameter
+ modifications.
+ cleanupParameters; Parameter values for pole cleanup. Available fields
+ are:
+ - 'boolCondition': if evaluation on pole returns False, pole is
+ removed;
+ - 'residueCheck': whether to apply residue check;
+ - 'residueNPoints': number of sample points to be used in residue
+ estimation;
+ - 'residueRadius': sample radius to be used in residue estimation
+ by Cauchy formula;
+ - 'residueTol': target tolerance to be used in residue estimation.
+ energyNormMatrix: Sparse matrix (in CSC format) assoociated to
+ w-weighted H10 norm.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ Q: Numpy 1D vector containing complex coefficients of approximant
+ denominator.
+ P: Numpy 2D vector whose columns are FE dofs of coefficients of
+ approximant numerator.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ """
+
+ extraApproxParameters = ["M", "N", "polyBasis"]
+ polyBasisParameters = ["CHEBYSHEV", "GAUSSLEGENDRE"]
+
+ def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
+ ks : "Numpy 1D array" = np.array([0]),
+ ws : "Numpy 1D array" = None, w : float = None,
+ approxParameters : dict = {}, plotSnapshots : bool = False,
+ verboseRobust : int = 1, cleanupParameters : dict = {}):
+ ROMApproximantLagrange.__init__(self, HFEngine, HSEngine, ks, w,
+ approxParameters, plotSnapshots)
+ if ws == None:
+ ws = np.ones(np.size(self.ks))
+ self.ws = ws
+ self.verboseRobust = verboseRobust
+ self.cleanupParameters = cleanupParameters
+
+ def parameterList(self) -> list:
+ """
+ List containing keys which are allowed in approxParameters.
+
+ Returns:
+ List of strings.
+ """
+ return (ROMApproximantLagrange.parameterList(self)
+ + ROMApproximantLagrangePade.extraApproxParameters)
+
+ @property
+ def k0(self):
+ """Dummy center of approximant (i.e. k0)."""
+ self.k0 = np.mean(self.ks)
+ return self._k0
+ @k0.setter
+ def k0(self, k0:bool):
+ self._k0 = k0
+
+ @property
+ def approxParameters(self):
+ """
+ Value of approximant parameters. Its assignment may change M, N and S.
+ """
+ return self._approxParameters
+ @approxParameters.setter
+ def approxParameters(self, approxParams):
+ if not hasattr(self, "approxParameters"):
+ self._approxParameters = {}
+ approxParameters = utilities.purgeDict(approxParams,
+ ROMApproximantLagrangePade.parameterList(self),
+ dictname = "approxParameters")
+ keyList = list(approxParameters.keys())
+ if "M" in keyList:
+ self.M = 0 #to avoid warnings if M and S are changed simultaneously
+ if "N" in keyList:
+ self.N = 0 #to avoid warnings if N and S are changed simultaneously
+ approxParametersCopy = utilities.purgeDict(approxParameters,
+ ROMApproximantLagrangePade.extraApproxParameters,
+ True, True)
+ ROMApproximantLagrange.approxParameters.fset(self,approxParametersCopy)
+ if "M" in keyList:
+ self.M = approxParameters["M"]
+ elif hasattr(self, "M"):
+ self.M = self.M
+ else:
+ self.M = 0
+ if "N" in keyList:
+ self.N = approxParameters["N"]
+ elif hasattr(self, "N"):
+ self.N = self.N
+ else:
+ self.N = 0
+ if "polyBasis" in keyList:
+ self.polyBasis = approxParameters["polyBasis"]
+ elif hasattr(self, "polyBasis"):
+ self.polyBasis = self.polyBasis
+ else:
+ self.polyBasis = "CHEBYSHEV"
+
+ @property
+ def M(self):
+ """Value of M. Its assignment may change S."""
+ return self._M
+ @M.setter
+ def M(self, M):
+ if M < 0: raise ArithmeticError("M must be non-negative.")
+ self._M = M
+ self._approxParameters["M"] = self.M
+ if hasattr(self, "S") and self.S < self.M + 1:
+ warnings.warn("Prescribed S is too small. Updating S to M + 1.")
+ self.S = self.M + 1
+
+ @property
+ def N(self):
+ """Value of N. Its assignment may change S."""
+ return self._N
+ @N.setter
+ def N(self, N):
+ if N < 0: raise ArithmeticError("N must be non-negative.")
+ self._N = N
+ self._approxParameters["N"] = self.N
+ if hasattr(self, "S") and self.S < self.N + 1:
+ warnings.warn("Prescribed S is too small. Updating S to N + 1.")
+ self.S = self.N + 1
+
+ @property
+ def S(self):
+ """Value of S."""
+ return self._S
+ @S.setter
+ def S(self, S):
+ if S <= 0: raise ArithmeticError("S must be positive.")
+ if hasattr(self, "S"): Sold = self.S
+ else: Sold = -1
+ vals, label = [0] * 2, {0:"M", 1:"N"}
+ if hasattr(self, "M"): vals[0] = self.M
+ if hasattr(self, "N"): vals[1] = self.N
+ idxmax = np.argmax(vals)
+ if vals[idxmax] + 1 > S:
+ warnings.warn("Prescribed S is too small. Updating S to {} + 1."\
+ .format(label[idxmax]))
+ self.Emax = vals[idxmax] + 1
+ else:
+ self._S = S
+ self._approxParameters["S"] = self.S
+ if Sold != self.S:
+ self.resetSamples()
+
+ @property
+ def polyBasis(self):
+ """Value of polyBasis."""
+ return self._polyBasis
+ @polyBasis.setter
+ def polyBasis(self, polyBasis):
+ if hasattr(self, "polyBasis"): polyBasisold = self.polyBasis
+ else: polyBasisold = -1
+ try:
+ polyBasis = polyBasis.upper().strip().replace(" ","")
+ if polyBasis not in self.polyBasisParameters: raise
+ except:
+ warnings.warn(("Prescribed polyBasis not recognized. Overriding to"
+ " 'CHEBYSHEV'."))
+ polyBasis = "CHEBYSHEV"
+ self._polyBasis = polyBasis
+ self._approxParameters["polyBasis"] = self.polyBasis
+ if polyBasisold != self.polyBasis:
+ self.resetSamples()
+
+ @property
+ def cleanupParameters(self):
+ """Value of cleanupParameters."""
+ return self._cleanupParameters
+ @cleanupParameters.setter
+ def cleanupParameters(self, cleanupParameters):
+ allowedCleanupKeys = ['boolCondition','residueCheck','residueNPoints',
+ 'residueRadius','residueTol']
+ cleanupKeys = cleanupParameters.keys()
+ if 'boolCondition' not in cleanupKeys:
+ cleanupParameters['boolCondition'] = lambda x : True
+ if 'residueCheck' not in cleanupKeys:
+ cleanupParameters['residueCheck'] = False
+ if 'residueNPoints' not in cleanupKeys:
+ cleanupParameters['residueNPoints'] = 2
+ if 'residueRadius' not in cleanupKeys:
+ cleanupParameters['residueRadius'] = 1e-5
+ if 'residueTol' not in cleanupKeys:
+ cleanupParameters['residueTol'] = 1e-4
+ cleanupParameters['boolCondition'] = np.vectorize(
+ cleanupParameters['boolCondition'])
+ self._cleanupParameters = {key : cleanupParameters[key]
+ for key in allowedCleanupKeys}
+
+ def setupApprox(self):
+ """
+ Compute Pade' interpolant.
+ SVD-based robust eigenvalue management.
+ """
+ S0 = self.S
+ M1 = self.M + 1
+ N1 = self.N + 1
+ if self.solSnapshots is None:
+ self.approxRadius = np.max(np.abs(self.k0 - self.ks))
+ if self.polyBasis == "CHEBYSHEV":
+ nodes, weights = np.polynomial.chebyshev.chebgauss(S0)
+ self.ws = weights * 2. / PI
+ elif self.polyBasis == "GAUSSLEGENDRE":
+ nodes, weights = np.polynomial.legendre.leggauss(S0)
+ self.ws = weights
+ self.ks = nodes * self.approxRadius + self.k0
+ self.ws = self.ws[:, None]
+ self.computeSnapshots()
+ if not self.checkComputedApprox():
+ Phi = np.zeros((S0, max(M1, N1)), dtype = np.complex)
+ Phi[:, 0] = np.ones((S0,)) / 2**.5
+ for j in range(1, max(M1, N1)):
+ c = np.zeros((j + 1,))
+ c[-1] = 1.
+ if self.polyBasis == "CHEBYSHEV":
+ Phi[:, j] = np.polynomial.chebyshev.chebval(
+ self.radiusPade(self.ks), c)
+ elif self.polyBasis == "GAUSSLEGENDRE":
+ Phi[:, j] = (j + .5) ** .5 * np.polynomial.legendre.legval(
+ self.radiusPade(self.ks), c)
+
+ if self.POD:
+ II = np.array(np.arange(0, S0**3, S0)
+ + np.kron(np.arange(S0), np.ones(S0)),
+ dtype = np.int)
+ Rtall = np.zeros(S0**3, dtype = np.complex)
+ Rtall[II] = np.reshape(self.RPOD.T, (S0**2,))
+ Rtall = np.reshape(Rtall, (S0, S0**2)).T
+
+ B = Rtall.dot(Phi[:, : N1])
+ Z = copy(B)
+ Z = np.reshape(Z.T, (S0 * (N1), S0)).T
+
+ for j in range(2):
+ Z = Z - Phi[:, : M1].dot(Phi[:, : M1].conj().T.dot(
+ np.multiply(self.ws, Z)))
+ Z = np.reshape(Z.T, (N1, S0**2)).T
+ else:
+ ker = self.solSnapshots.conj().T.dot(self.energyNormMat.dot(
+ self.solSnapshots))
+ WPhi = np.reshape(np.multiply(self.ws, Phi[:, : M1]).T,
+ (S0 * M1)).conj()[:, None]
+ Y = np.multiply(WPhi, np.kron(np.ones((M1, 1)), Phi[:, : N1]))
+ Ylarge = np.reshape(Y.T, (M1 * N1, S0)).T
+
+ B = np.reshape(ker.dot(Ylarge).T, (N1, M1 * S0)).T
+ D = np.multiply(self.ws, np.diag(ker)[:, None])
+ D = Phi[:, : N1].conj().T.dot(np.multiply(D, Phi[:, : N1]))
+
+ Z = B.conj().T.dot(Y) - D
+ _, ev, eV = np.linalg.svd(Z, full_matrices=False)
+ eV = eV.conj().T
+ phi = eV[:, np.argmin(ev)]
+
+ if self.POD:
+ c = Phi[:, : M1].conj().T.dot(np.multiply(self.ws,
+ np.reshape(B.dot(phi).T, (S0, S0)).T))
+ else:
+ c = np.reshape(Y.dot(phi).T, (M1, S0))
+
+ polybase = np.zeros((max(M1, N1), max(M1, N1)))
+ polybase[0, 0] = 1
+ polybase[1, 1] = 1
+ for j in range(2, max(M1, N1)):
+ if self.polyBasis == "CHEBYSHEV":
+ polybase[1 : j + 1, j] = 2 * polybase[: j, j - 1]
+ polybase[: j + 1, j] = (polybase[: j + 1, j]
+ - polybase[: j + 1, j - 2])
+ elif self.polyBasis == "GAUSSLEGENDRE":
+ polybase[1 : j + 1, j] = ((2 * j - 1) / j
+ * polybase[: j, j - 1])
+ polybase[: j + 1, j] = (polybase[: j + 1, j]
+ - ((j - 1) / j * polybase[: j + 1, j - 2]))
+ if self.polyBasis == "CHEBYSHEV":
+ polybase[0, 0] = .5 ** .5
+ elif self.polyBasis == "GAUSSLEGENDRE":
+ polybase = np.multiply(np.power(
+ np.arange(.5, max(M1, N1)), .5), polybase)
+ self.P = self.solSnapshots.dot(polybase[: M1, : M1].dot(c).T)
+ self.Q = polybase[: N1, : N1].dot(phi)
+
+ self.approxParameters = {"N" : self.N, "M" : self.M, "S" : S0}
+ self.lastApproxParameters = copy(self.approxParameters)
+
+ def _cleanup(self):
+ """Cleanup Pade' denominator by removing unwanted poles."""
+ if self.N == 0: return
+ poles = np.roots(self.Q[::-1]) + self.k0
+ NpolesOld = len(poles)
+ poles = poles[self.cleanupParameters['boolCondition'](poles)]
+
+ if self.cleanupParameters['residueCheck']:
+ resR = self.cleanupParameters['residueRadius']
+ resTol = self.cleanupParameters['residueTol']
+ residues = np.zeros_like(poles)
+ for l, pole in enumerate(poles):
+ resV = np.zeros(self.solDerivatives.shape[0],
+ dtype = np.complex)
+ NPoints = self.cleanupParameters['residueNPoints']
+ for theta in 2 * PI * np.arange(NPoints) / NPoints:
+ deltapole = resR * np.exp(1.j * theta)
+ ksample = pole + deltapole
+ self.solveHF(ksample)
+ resV = deltapole / NPoints * (resV + self.uHF)
+ residues[l] = np.abs(resV.dot(
+ self.energyNormMatrix.dot(resV).conj())) ** .5
+ poles = poles[residues >= resTol]
+
+ NpolesNew = len(poles)
+ if NpolesOld > NpolesNew:
+ if self.verboseRobust >= 1:
+ sSing = "s" * (NpolesOld - NpolesNew > 1)
+ print(("Identified {0} pole{1} to be removed.\n"
+ "Reducing N from {2} to {3}.").format(
+ NpolesOld - NpolesNew, sSing, NpolesOld, NpolesNew))
+ self.N = NpolesNew
+ newQ = np.polyfit(np.append(poles - self.k0, 0.),
+ np.append(np.zeros(NpolesNew), 1.), NpolesNew)
+ self.Q = newQ[::-1] / np.linalg.norm(newQ)
+
+ def radiusPade(self, k:complex):
+ """
+ Compute translated radius to be plugged into Pade' approximant.
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ Translated radius to be plugged into Pade' approximant.
+ """
+ return (k - self.k0) / self.approxRadius
+
+ def evalApprox(self, k:complex) -> ("Fenics function", "Fenics function"):
+ """
+ Evaluate Pade' approximant at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ Real and imaginary parts of approximant.
+ """
+ self.setupApprox()
+ powerlist = np.power(self.radiusPade(k), range(max(self.M, self.N) + 1))
+ self.uApp = (self.P.dot(powerlist[:self.M + 1])
+ / self.Q.dot(powerlist[:self.N + 1]))
+ return self.uApp
+
+ def getPoles(self, centered : bool = False) -> "numpy 1D array":
+ """
+ Obtain approximant poles.
+
+ Args:
+ centered(optional): Whether to return pole positions relative to
+ approximation center. Defaults to False.
+
+ Returns:
+ Numpy complex vector of poles.
+ """
+ self.setupApprox()
+ return np.roots(self.Q[::-1]) * self.approxRadius + self.k0 * centered
+
+
\ No newline at end of file
diff --git a/main/ROMApproximantLagrangeRB.py b/main/ROMApproximantLagrangeRB.py
new file mode 100644
index 0000000..47aab74
--- /dev/null
+++ b/main/ROMApproximantLagrangeRB.py
@@ -0,0 +1,301 @@
+#!/usr/bin/python
+
+from copy import copy
+import warnings
+import numpy as np
+import scipy as sp
+import utilities
+from ROMApproximantLagrange import ROMApproximantLagrange
+
+PI = np.pi
+
+class ROMApproximantLagrangeRB(ROMApproximantLagrange):
+ """
+ ROM RB approximant computation for parametric problems.
+
+ Args:
+ HFEngine: HF problem solver. Should have members:
+ - energyNormMatrix: assemble sparse matrix (in CSC format)
+ associated to weighted H10 norm;
+ - problemData: list of HF problem data (excluding k);
+ - setProblemData: set HF problem data and k to prescribed values;
+ - solve: get HF solution as complex numpy vector.
+ HSEngine: Hilbert space general purpose engine. Should have members:
+ - norm: compute Hilbert norm of expression represented as complex
+ numpy vector;
+ - plot: plot Hilbert expression represented as complex numpy vector.
+ ks(optional): Array of snapshot parameters. Defaults to np.array([0]).
+ ws(optional): Array of snapshots weights. Defaults to uniform = 1.
+ w(optional): Weight for norm computation. If None, set to
+ Re(np.mean(ks)). Defaults to None.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True;
+ - 'S': total number of samples current approximant relies upon;
+ defaults to 2;
+ - 'R': rank for Galerkin projection; defaults to S;
+ - 'nodesType': sampling node type; available values are:
+ - 'CHEBYSHEV': Chebyshev nodes;
+ - 'GAUSSLEGENDRE': Gauss-Legendre nodes;
+ - 'MANUAL': manual selection;
+ defaults to 'MANUAL'.
+ Defaults to empty dict.
+ plotSnapshots(optional): Whether to plot snapshots of the Helmholtz
+ solution map. Defaults to False.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ HSEngine: Hilbert space general purpose engine.
+ solSnapshots: Array whose columns are FE dofs of snapshots of solution
+ map.
+ k0: Default parameter.
+ ks: Array of snapshot parameters.
+ ws: Array of snapshots weights.
+ w: Weight for norm computation.
+ approxRadius: Dummy radius of approximant (i.e. distance from k0 to
+ farthest sample point).
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots;
+ - 'S': total number of samples current approximant relies upon;
+ - 'R': rank for Galerkin projection;
+ - 'nodesType': sampling node type.
+ S: Number of solution snapshots over which current approximant is
+ based upon.
+ R: Rank for Galerkin projection.
+ nodesType: Sampling node type.
+ plotSnapshots: Whether to plot snapshots of the Helmholtz solution map.
+ POD: Whether to compute POD of snapshots.
+ energyNormMatrix: Sparse matrix (in CSC format) assoociated to
+ w-weighted H10 norm.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ """
+
+ extraApproxParameters = ["R", "nodesType"]
+ nodesTypeParameters = ["CHEBYSHEV", "GAUSSLEGENDRE", "MANUAL"]
+
+ def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
+ ks : "Numpy 1D array" = np.array([0]),
+ ws : "Numpy 1D array" = None, w : float = None,
+ approxParameters : dict = {}, plotSnapshots : bool = False):
+ ROMApproximantLagrange.__init__(self, HFEngine, HSEngine, ks, w,
+ approxParameters, plotSnapshots)
+ if ws == None:
+ ws = np.ones(np.size(self.ks))
+ self.ws = ws
+ self.solLifting = self.HFEngine.liftDirichletData()
+
+ @property
+ def k0(self):
+ """Dummy center of approximant (i.e. k0)."""
+ self.k0 = np.mean(self.ks)
+ return self._k0
+ @k0.setter
+ def k0(self, k0:bool):
+ self._k0 = k0
+
+ def resetSamples(self):
+ """Reset samples."""
+ ROMApproximantLagrange.resetSamples(self)
+ self.projMat = None
+
+ def parameterList(self) -> list:
+ """
+ List containing keys which are allowed in approxParameters.
+
+ Returns:
+ List of strings.
+ """
+ return (ROMApproximantLagrange.parameterList(self)
+ + ROMApproximantLagrangeRB.extraApproxParameters)
+
+ @property
+ def approxParameters(self):
+ """
+ Value of approximant parameters. Its assignment may change M, N and S.
+ """
+ return self._approxParameters
+ @approxParameters.setter
+ def approxParameters(self, approxParams):
+ if not hasattr(self, "approxParameters"):
+ self._approxParameters = {}
+ approxParameters = utilities.purgeDict(approxParams,
+ ROMApproximantLagrangeRB.parameterList(self),
+ dictname = "approxParameters")
+ keyList = list(approxParameters.keys())
+ approxParametersCopy = utilities.purgeDict(approxParameters,
+ ROMApproximantLagrangeRB.extraApproxParameters,
+ True, True)
+ ROMApproximantLagrange.approxParameters.fset(self,approxParametersCopy)
+ if "R" in keyList:
+ self.R = approxParameters["R"]
+ elif hasattr(self, "R"):
+ self.R = self.R
+ else:
+ self.R = self.S
+ if "nodesType" in keyList:
+ self.nodesType = approxParameters["nodesType"]
+ elif hasattr(self, "nodesType"):
+ self.nodesType = self.nodesType
+ else:
+ self.nodesType = "MANUAL"
+
+ @property
+ def R(self):
+ """Value of R. Its assignment may change S."""
+ return self._R
+ @R.setter
+ def R(self, R):
+ if R < 0: raise ArithmeticError("R must be non-negative.")
+ self._R = R
+ self._approxParameters["R"] = self.R
+ if hasattr(self, "S") and self.S < self.R:
+ warnings.warn("Prescribed S is too small. Updating S to R.")
+ self.S = self.R
+
+ @property
+ def nodesType(self):
+ """Value of nodesType."""
+ return self._nodesType
+ @nodesType.setter
+ def nodesType(self, nodesType):
+ if hasattr(self, "nodesType"): nodesTypeold = self.nodesType
+ else: nodesTypeold = -1
+ try:
+ nodesType = nodesType.upper().strip().replace(" ","")
+ if nodesType not in self.nodesTypeParameters: raise
+ except:
+ warnings.warn(("Prescribed nodesType not recognized. Overriding to"
+ " 'MANUAL'."))
+ nodesType = "MANUAL"
+ self._nodesType = nodesType
+ self._approxParameters["nodesType"] = self.nodesType
+ if nodesTypeold != self.nodesType:
+ self.resetSamples()
+
+ def manageSnapshots(self, u:"2-tuple of Fenics function", pos:int):
+ """
+ Post-process snapshots of solution map.
+
+ Any specialization should include something like
+ self.solSnapshots[:, pos] = (np.array(u[0].vector()[:])
+ + 1.j * np.array(u[1].vector()[:]))
+
+ Args:
+ u: 2-tuple containing real and imaginary parts of FE dofs of
+ snapshot.
+ pos: Derivative index.
+ """
+ if pos == 0:
+ self.As, self.bs = self.HFEngine.getLSBlocks()
+ u = u - self.solLifting
+ return ROMApproximantLagrange.manageSnapshots(self, u, pos)
+
+ def setupApprox(self):
+ """
+ Compute RB projection matrix.
+ See ``Householder triangularization of a quasimatrix'', L.Trefethen,
+ 2008 for QR algorithm.
+ """
+ if self.solSnapshots is None:
+ if self.nodesType == "MANUAL" and len(self.ks) != self.S:
+ warnings.warn(("Number of prescribed nodes different from S. "
+ "Overriding S to len(ks)"))
+ self.S = len(self.ks)
+ self.approxRadius = np.max(np.abs(self.k0 - self.ks))
+ if self.nodesType != "MANUAL":
+ if self.nodesType == "CHEBYSHEV":
+ nodes, weights = np.polynomial.chebyshev.chebgauss(self.S)
+ self.ws = weights * 2. / PI
+ elif self.nodesType == "GAUSSLEGENDRE":
+ nodes, weights = np.polynomial.legendre.leggauss(self.S)
+ self.ws = weights
+ self.ks = nodes * self.approxRadius + self.k0
+ self.ws = self.ws[:, None]
+ self.computeSnapshots()
+ if not self.checkComputedApprox():
+ if self.POD:
+ U, _, _ = np.linalg.svd(self.RPOD, full_matrices = False)
+ self.projMat = self.solSnapshots.dot(U[:, : self.R])
+ else:
+ self.projMat = self.solSnapshots[:, : self.R]
+ self.assembleReducedSystem()
+ self.lastApproxParameters = copy(self.approxParameters)
+
+ def assembleReducedSystem(self):
+ """Build affine blocks of RB linear system through projections."""
+ projMatH = self.projMat.T.conjugate()
+ self.ARBs = [None] * len(self.As)
+ self.bRBs = [None] * max(len(self.As), len(self.bs))
+ for j in range(len(self.As)):
+ self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat))
+ if j < len(self.bs):
+ self.bRBs[j] = projMatH.dot(self.bs[j]
+ - self.As[j].dot(self.solLifting))
+ else:
+ self.bRBs[j] = - projMatH.dot(self.As[j].dot(self.solLifting))
+ for j in range(len(self.As), len(self.bs)):
+ self.bRBs[j] = projMatH.dot(self.bs[j])
+
+ def solveReducedSystem(self, k:complex) -> "Numpy 1D array":
+ """
+ Solve RB linear system.
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ Solution of RB linear system.
+ """
+ self.setupApprox()
+ ARBk = self.ARBs[0][:self.R,:self.R]
+ bRBk = self.bRBs[0][:self.R]
+ for j in range(1, len(self.ARBs)):
+ ARBk = ARBk + np.power(k, j) * self.ARBs[j][:self.R, :self.R]
+ for j in range(1, len(self.bRBs)):
+ bRBk = bRBk + np.power(k, j) * self.bRBs[j][:self.R]
+ return self.projMat[:, :self.R].dot(np.linalg.solve(ARBk, bRBk))
+
+ def evalApprox(self, k:complex) -> ("Fenics function", "Fenics function"):
+ """
+ Evaluate RB approximant at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ Real and imaginary parts of approximant.
+ """
+ self.setupApprox()
+ self.uApp = self.solLifting + self.solveReducedSystem(k)
+ return self.uApp
+
+ def getPoles(self, centered : bool = False) -> "numpy 1D array":
+ """
+ Obtain approximant poles.
+
+ Args:
+ centered(optional): Whether to return pole positions relative to
+ approximation center. Defaults to False.
+
+ Returns:
+ Numpy complex vector of poles.
+ """
+ self.setupApprox()
+ A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1),
+ dtype = np.complex)
+ B = np.zeros_like(A)
+ A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0]
+ for j in range(len(self.ARBs) - 1):
+ Aj = self.ARBs[j + 1]
+ B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj
+ II = np.arange(self.ARBs[0].shape[0],
+ self.ARBs[0].shape[0] * (len(self.ARBs) - 1))
+ B[II, II - self.ARBs[0].shape[0]] = 1.
+ return sp.linalg.eigvals(A, B) - self.k0 * (not centered)
diff --git a/main/ROMApproximantSweeper.py b/main/ROMApproximantSweeper.py
new file mode 100644
index 0000000..0378b1d
--- /dev/null
+++ b/main/ROMApproximantSweeper.py
@@ -0,0 +1,257 @@
+#!/usr/bin/python
+
+import os
+import csv
+import warnings
+import numpy as np
+from context import utilities
+
+class ROMApproximantSweeper:
+ """
+ ROM approximant sweeper.
+
+ Args:
+ ROMEngine(optional): ROMApproximant class. Defaults to None.
+ ktars(optional): Array of parameter values to sweep. Defaults to empty
+ array.
+ params(optional): List of parameter settings (each as a dict) to
+ explore. Defaults to single empty set.
+ mostExpensive(optional): String containing label of most expensive
+ step, to be executed fewer times. Allowed options are 'HF' and
+ 'Approx'. Defaults to 'HF'.
+
+ Attributes:
+ ROMEngine: ROMApproximant class.
+ ktars: Array of parameter values to sweep.
+ params: List of parameter settings (each as a dict) to explore.
+ mostExpensive: String containing label of most expensive step, to be
+ executed fewer times.
+ """
+
+ allowedOutputs = ["HFNorm", "AppNorm", "AppError"]
+
+ def __init__(self, ROMEngine : 'ROMApproximant' = None,
+ ktars : 'numpy 1D array' = np.array([]),
+ params : "List of dicts" = [{}],
+ mostExpensive : str = "HF"):
+ self.ROMEngine = ROMEngine
+ self.ktars = ktars
+ self.params = params
+ self.mostExpensive = mostExpensive
+
+ def name(self) -> str:
+ """Approximant label."""
+ return self.__class__.__name__
+
+ @property
+ def mostExpensive(self):
+ """Value of mostExpensive."""
+ return self._mostExpensive
+ @mostExpensive.setter
+ def mostExpensive(self, mostExpensive:str):
+ mostExpensive = mostExpensive.upper()
+ if mostExpensive not in ["HF", "APPROX"]:
+ warnings.warn(("Value of mostExpensive not recognized. Overriding "
+ "to 'HF'."))
+ mostExpensive = "HF"
+ self._mostExpensive = mostExpensive
+
+ def checkValues(self) -> bool:
+ """Check if sweep can be performed."""
+ if self.ROMEngine is None:
+ warnings.warn("ROMEngine is missing. Aborting.")
+ return False
+ if len(self.ktars) == 0:
+ warnings.warn("Empty target parameter vector. Aborting.")
+ return False
+ if len(self.params) == 0:
+ warnings.warn("Empty method parameters vector. Aborting.")
+ return False
+ return True
+
+ def sweep(self, filename : str = "./out", outputs : list = [],
+ verbose : int = 1):
+ if not self.checkValues(): return
+ try:
+ if outputs.upper() == "ALL":
+ outputs = self.allowedOutputs + ["poles"]
+ except:
+ if len(outputs) == 0:
+ outputs = self.allowedOutputs
+ outputs = utilities.purgeList(outputs,
+ self.allowedOutputs + ["poles"],
+ listname = "outputs list")
+ poles = ("poles" in outputs)
+ if len(outputs) == 0:
+ warnings.warn("Empty outputs. Aborting.")
+ return
+ outParList = self.ROMEngine.parameterList()
+ Nparams = len(self.params)
+ allowedParams = self.ROMEngine.parameterList()
+
+ while os.path.exists(filename):
+ filename = filename + "{}".format(np.random.randint(10))
+ append_write = "w"
+ initial_row = (outParList + ["kRe", "kIm"]
+ + [x for x in self.allowedOutputs if x in outputs]
+ + ["type"] + ["poles"] * poles)
+ with open(filename, append_write, buffering = 1) as fout:
+ writer = csv.writer(fout, delimiter=",")
+ writer.writerow(initial_row)
+
+ if self.mostExpensive == "HF":
+ outerSet = self.ktars
+ innerSet = self.params
+ elif self.mostExpensive == "APPROX":
+ outerSet = self.params
+ innerSet = self.ktars
+
+ for outerIdx, outerPar in enumerate(outerSet):
+ if self.mostExpensive == "HF":
+ i, ktar = outerIdx, outerPar
+ elif self.mostExpensive == "APPROX":
+ j, par = outerIdx, outerPar
+ self.ROMEngine.approxParameters = {k: par[k] for k in\
+ par.keys() & allowedParams}
+ self.ROMEngine.setupApprox()
+
+ for innerIdx, innerPar in enumerate(innerSet):
+ if self.mostExpensive == "APPROX":
+ i, ktar = innerIdx, innerPar
+ elif self.mostExpensive == "HF":
+ j, par = innerIdx, innerPar
+ self.ROMEngine.approxParameters = {k: par[k] for k in\
+ par.keys() & allowedParams}
+ self.ROMEngine.setupApprox()
+
+ if verbose >= 1:
+ print("Set {}/{}\tk_{} = {:.10f}{}"\
+ .format(j+1, Nparams, i, ktar, " " * 15),
+ end="\r")
+
+ outData = []
+ if "HFNorm" in outputs:
+ outData = outData + [self.ROMEngine.HFNorm(ktar)]
+ if "AppNorm" in outputs:
+ outData = outData + [self.ROMEngine.approxNorm(ktar)]
+ if "AppError" in outputs:
+ outData = outData + [self.ROMEngine.approxError(ktar)]
+ writeData = []
+ for parn in outParList:
+ writeData = (writeData
+ + [self.ROMEngine.approxParameters[parn]])
+ writeData = (writeData + [ktar.real, ktar.imag]
+ + outData + [self.ROMEngine.name()])
+ if poles:
+ writeData = writeData + list(self.ROMEngine.getPoles())
+ writer.writerow(str(x) for x in writeData)
+
+ if verbose >= 1:
+ if self.mostExpensive == "APPROX":
+ print("Set {}/{}\tdone{}".format(j+1, Nparams, " "*25))
+ elif self.mostExpensive == "HF":
+ print("Point k_{} = {:.10f}\tdone{}".format(i, ktar,
+ " " * 25))
+ self.filename = filename
+ return self.filename
+
+ def read(self, filename:str, restrictions : dict = {},
+ outputs : list = []) -> dict:
+ """
+ Execute a query on a custom format CSV.
+
+ Args:
+ filename: CSV filename.
+ restrictions(optional): Parameter configurations to output.
+ Defaults to empty dictionary, i.e. output all.
+ outputs(optional): Values to output. Defaults to empty list, i.e.
+ no output.
+
+ Returns:
+ Dictionary of desired results, with a key for each entry of
+ outputs, and a numpy 1D array as corresponding value.
+ """
+ def pairKFromCSV(filename:str, kTarget:complex) -> "2-tuple of str":
+ """
+ Find complex point in CSV closer to a prescribed value.
+
+ Args:
+ filename: CSV filename.
+ zTarget: Target complex value.
+
+ Returns:
+ Strings containing real and imaginary part of desired value, in
+ the same format as in the CSV file.
+ """
+ ktarsF = np.array([], dtype = complex)
+ kRetarsF = np.array([], dtype = complex)
+ kImtarsF = np.array([], dtype = complex)
+ with open(filename, 'r') as f:
+ reader = csv.reader(f, delimiter=',')
+ header = next(reader)
+ kReindex = header.index('kRe')
+ kImindex = header.index('kIm')
+ for row in reader:
+ try:
+ if row[kReindex] not in [" ", ""]:
+ kRetarsF = np.append(kRetarsF, row[kReindex])
+ kImtarsF = np.append(kImtarsF, row[kImindex])
+ ktarsF = np.append(ktarsF, float(row[kReindex])
+ + 1.j * float(row[kImindex]))
+ except:
+ pass
+ optimalIndex = np.argmin(np.abs(ktarsF - kTarget))
+ return [kRetarsF[optimalIndex], kImtarsF[optimalIndex]]
+
+ with open(filename, 'r') as f:
+ reader = csv.reader(f, delimiter=',')
+ header = next(reader)
+ restrIndices, outputIndices, outputData = {}, {}, {}
+ for key in restrictions.keys():
+ try:
+ restrIndices[key] = header.index(key)
+ if not isinstance(restrictions[key], list):
+ restrictions[key] = [restrictions[key]]
+ restrictions[key] = [str(x) for x in restrictions[key]]
+ except:
+ warnings.warn("Ignoring key {} from restrictions".format(
+ key))
+
+ if 'kRe' in restrIndices.keys() or 'kIm' in restrIndices.keys():
+ if 'kRe' not in restrIndices.keys():
+ restrIndices['kRe'] = header.index('kRe')
+ restrictions['kRe'] = [0.] * len(restrictions['kIm'])
+ elif 'kIm' not in restrIndices.keys():
+ restrIndices['kIm'] = header.index('kIm')
+ restrictions['kIm'] = [0.] * len(restrictions['kRe'])
+ elif len(restrictions['kRe']) != len(restrictions['kIm']):
+ raise Exception(("The lists of values for kRe and kIm "
+ "must have the same length."))
+ for i in range(len(restrictions['kRe'])):
+ k = (1.0 * float(restrictions['kRe'][i])
+ + 1.j * float(restrictions['kIm'][i]))
+ restrictions['kRe'][i], restrictions['kIm'][i] =\
+ pairKFromCSV(filename, k)
+ for key in outputs:
+ try:
+ outputIndices[key] = header.index(key)
+ outputData[key] = np.array([])
+ except:
+ warnings.warn("Ignoring key {} from outputs".format(key))
+
+ for row in reader:
+ if all([row[restrIndices[key]] in restrictions[key]\
+ for key in restrictions.keys()]):
+ for key in outputIndices.keys():
+ try:
+ val = row[outputIndices[key]]
+ val = int(val)
+ except:
+ try:
+ val = float(val)
+ except:
+ val = np.nan
+ finally:
+ outputData[key] = np.append(outputData[key], val)
+ return outputData
+
\ No newline at end of file
diff --git a/main/ROMApproximantTaylor.py b/main/ROMApproximantTaylor.py
new file mode 100644
index 0000000..f554993
--- /dev/null
+++ b/main/ROMApproximantTaylor.py
@@ -0,0 +1,264 @@
+#!/usr/bin/python
+
+import warnings
+import numpy as np
+import utilities
+from ROMApproximant import ROMApproximant
+
+
+class ROMApproximantTaylor(ROMApproximant):
+ """
+ ROM single-point approximant computation for parametric problems.
+
+ Args:
+ HFEngine: HF problem solver. Should have members:
+ - energyNormMatrix: assemble sparse matrix (in CSC format)
+ associated to weighted H10 norm;
+ - problemData: list of HF problem data (excluding k);
+ - setProblemData: set HF problem data and k to prescribed values;
+ - setupDerivativeComputation: setup HF problem data to compute j-th
+ solution derivative at k0;
+ - solve: get HF solution as complex numpy vector.
+ HSEngine: Hilbert space general purpose engine. Should have members:
+ - norm: compute Hilbert norm of expression represented as complex
+ numpy vector;
+ - plot: plot Hilbert expression represented as complex numpy vector.
+ k0(optional): Default parameter. Defaults to 0.
+ w(optional): Weight for norm computation. If None, set to Re(k0).
+ Defaults to None.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True;
+ - 'E': total number of derivatives current approximant relies upon;
+ defaults to Emax;
+ - 'Emax': total number of derivatives of solution map to be
+ computed; defaults to E;
+ - 'sampleType': label of sampling type; available values are:
+ - 'ARNOLDI': orthogonalization of solution derivatives through
+ Arnoldi algorithm;
+ - 'KRYLOV': standard computation of solution derivatives.
+ Defaults to 'KRYLOV'.
+ Defaults to empty dict.
+ plotDer(optional): Whether to plot derivatives of the Helmholtz
+ solution map. Defaults to False.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ HSEngine: Hilbert space general purpose engine.
+ solDerivatives: Array whose columns are FE dofs of solution
+ derivatives.
+ k0: Default parameter.
+ w: Weight for norm computation.
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots;
+ - 'E': total number of derivatives current approximant relies upon;
+ - 'Emax': total number of derivatives of solution map to be
+ computed;
+ - 'sampleType': label of sampling type.
+ E: Number of solution derivatives over which current approximant is
+ based upon.
+ Emax: Total number of solution derivatives to be computed.
+ sampleType: Label of sampling type.
+ plotDer: Whether to plot derivatives of the Helmholtz solution map.
+ initialHFData: HF problem initial data.
+ energyNormMatrix: Sparse matrix (in CSC format) assoociated to
+ w-weighted H10 norm.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ """
+
+ extraApproxParameters = ["E", "Emax", "sampleType"]
+
+ def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
+ k0 : complex = 0, w : float = None,
+ approxParameters : dict = {}, plotDer : bool = False):
+ self.plotDer = plotDer
+ ROMApproximant.__init__(self, HFEngine, HSEngine,
+ k0, w, approxParameters)
+
+ def parameterList(self) -> list:
+ """
+ List containing keys which are allowed in approxParameters.
+
+ Returns:
+ List of strings.
+ """
+ return (ROMApproximant.parameterList(self)
+ + ROMApproximantTaylor.extraApproxParameters)
+
+ @property
+ def approxParameters(self):
+ """
+ Value of approximant parameters. Its assignment may change E and Emax.
+ """
+ return self._approxParameters
+ @approxParameters.setter
+ def approxParameters(self, approxParams):
+ if not hasattr(self, "approxParameters"):
+ self._approxParameters = {}
+ approxParameters = utilities.purgeDict(approxParams,
+ ROMApproximantTaylor.parameterList(self),
+ dictname = "approxParameters")
+ keyList = list(approxParameters.keys())
+ approxParametersCopy = utilities.purgeDict(approxParameters,
+ ROMApproximantTaylor.extraApproxParameters,
+ True, True)
+ ROMApproximant.approxParameters.fset(self, approxParametersCopy)
+ if "E" in keyList:
+ self._E = approxParameters["E"]
+ self._approxParameters["E"] = self.E
+ if "Emax" in keyList:
+ self.Emax = approxParameters["Emax"]
+ else:
+ if not hasattr(self, "Emax"):
+ self.Emax = self.E
+ else:
+ self.Emax = self.Emax
+ else:
+ if "Emax" in keyList:
+ self._E = approxParameters["Emax"]
+ self._approxParameters["E"] = self.E
+ self.Emax = self.E
+ else:
+ if not (hasattr(self, "Emax") and hasattr(self, "E")):
+ raise Exception("At least one of E and Emax must be set.")
+ if "sampleType" in keyList:
+ self.sampleType = approxParameters["sampleType"]
+ elif hasattr(self, "sampleType"):
+ self.sampleType = self.sampleType
+ else:
+ self.sampleType = "KRYLOV"
+
+ @property
+ def E(self):
+ """Value of E. Its assignment may change Emax."""
+ return self._E
+ @E.setter
+ def E(self, E):
+ if E < 0: raise ArithmeticError("E must be non-negative.")
+ self._E = E
+ self._approxParameters["E"] = self.E
+ if hasattr(self, "Emax") and self.Emax < self.E:
+ warnings.warn("Prescribed E is too large. Updating Emax to E.")
+ self.Emax = self.E
+
+ @property
+ def Emax(self):
+ """Value of Emax. Its assignment may reset computed derivatives."""
+ return self._Emax
+ @Emax.setter
+ def Emax(self, Emax):
+ if Emax < 0: raise ArithmeticError("Emax must be non-negative.")
+ if hasattr(self, "Emax"): EmaxOld = self.Emax
+ else: EmaxOld = -1
+ self._Emax = Emax
+ if hasattr(self, "E") and self.Emax < self.E:
+ warnings.warn("Prescribed Emax is too small. Updating Emax to E.")
+ self.Emax = self.E
+ else:
+ self._approxParameters["Emax"] = self.Emax
+ if EmaxOld >= self.Emax and self.solDerivatives is not None:
+ self.solDerivatives = self.solDerivatives[:, : self.Emax + 1]
+ if hasattr(self, "HArnoldi"):
+ self.HArnoldi = self.HArnoldi[: self.Emax + 1,
+ : self.Emax + 1]
+ self.RArnoldi = self.RArnoldi[: self.Emax + 1,
+ : self.Emax + 1]
+ else:
+ self.resetSamples()
+
+ @property
+ def sampleType(self):
+ """Value of sampleType."""
+ return self._sampleType
+ @sampleType.setter
+ def sampleType(self, sampleType):
+ if hasattr(self, "sampleType"): sampleTypeOld = self.sampleType
+ else: sampleTypeOld = -1
+ try:
+ sampleType = sampleType.upper().strip().replace(" ","")
+ if sampleType not in ["ARNOLDI", "KRYLOV"]: raise
+ if sampleType == "ARNOLDI" and not self.POD:
+ warnings.warn(("Prescribed sampleType not compatible with POD "
+ "value. Overriding to 'KRYLOV'."))
+ sampleType = "KRYLOV"
+ self._sampleType = sampleType
+ except:
+ warnings.warn(("Prescribed sampleType not recognized. Overriding "
+ "to 'KRYLOV'."))
+ self.sampleType = "KRYLOV"
+ self._approxParameters["sampleType"] = self.sampleType
+ if sampleTypeOld != self.sampleType:
+ self.resetSamples()
+
+ def resetSamples(self):
+ """Reset samples."""
+ self.solDerivatives = None
+ if hasattr(self, "HArnoldi"):
+ del self.HArnoldi
+ if hasattr(self, "RArnoldi"):
+ del self.RArnoldi
+
+ def computeDerivatives(self):
+ """
+ Compute derivatives of solution map starting from order 0.
+ """
+ if self.solDerivatives is None:
+ up = np.zeros(self.HSEngine.V.dim())
+ upp = up
+ for j in range(self.Emax + 1):
+ if j == 0:
+ self.HFEngine.setProblemData(self.initialHFData, self.k0)
+ else:
+ self.HFEngine.setupDerivativeComputation(j, up, upp)
+ uj = self.HFEngine.solve(True)
+ if self.plotDer:
+ self.HSEngine.plot(uj, name = "u_{0}".format(j))
+ upp = up
+ up = self.manageDerivatives(uj, j)
+
+ def manageDerivatives(self, u:"numpy 1D array",
+ pos:int) -> ("numpy 1D array"):
+ """
+ Store derivatives of solution map.
+
+ Args:
+ u: solution derivative as numpy complex vector;
+ pos: Derivative index.
+
+ Returns:
+ Adjusted solution derivative as numpy complex vector.
+ """
+ if pos == 0:
+ self.solDerivatives = np.empty((u.shape[0], self.Emax + 1),
+ dtype = np.complex)
+ self.solDerivatives[:, pos] = u
+ if self.sampleType == "ARNOLDI":
+ if pos == 0:
+ self.HArnoldi = np.zeros((self.Emax + 1, self.Emax + 1),
+ dtype = np.complex)
+ self.RArnoldi = np.zeros((self.Emax + 1, self.Emax + 1),
+ dtype = np.complex)
+ self.HArnoldi[: pos, pos] = self.solDerivatives[:, pos].conj().dot(
+ self.energyNormMatrix.dot(self.solDerivatives[:, : pos])
+ ).conj()
+ self.solDerivatives[:, pos] = (self.solDerivatives[:, pos]
+ - self.solDerivatives[:, : pos].dot(
+ self.HArnoldi[: pos, pos]))
+ self.HArnoldi[pos, pos] = (self.solDerivatives[:, pos].conj().dot(
+ self.energyNormMatrix.dot(self.solDerivatives[:, pos])))**.5
+ self.solDerivatives[:, pos] = (self.solDerivatives[:, pos]
+ / self.HArnoldi[pos, pos])
+ if pos == 0:
+ self.RArnoldi[pos, pos] = self.HArnoldi[pos, pos]
+ else:
+ self.RArnoldi[:pos+1, pos] = self.HArnoldi[: pos+1, 1 : pos+1]\
+ .dot(self.RArnoldi[: pos, pos - 1])
+ return self.solDerivatives[:, pos]
+
diff --git a/main/ROMApproximantTaylorPade.py b/main/ROMApproximantTaylorPade.py
new file mode 100644
index 0000000..e5b6d87
--- /dev/null
+++ b/main/ROMApproximantTaylorPade.py
@@ -0,0 +1,552 @@
+#!/usr/bin/python
+
+from __future__ import print_function
+from copy import copy
+import warnings
+import numpy as np
+import utilities
+from ROMApproximantTaylor import ROMApproximantTaylor
+
+PI = np.pi
+
+class ROMApproximantTaylorPade(ROMApproximantTaylor):
+ """
+ ROM single-point fast Pade' approximant computation for parametric
+ problems.
+
+ Args:
+ HFEngine: HF problem solver. Should have members:
+ - energyNormMatrix: assemble sparse matrix (in CSC format)
+ associated to weighted H10 norm;
+ - problemData: list of HF problem data (excluding k);
+ - setProblemData: set HF problem data and k to prescribed values;
+ - setupDerivativeComputation: setup HF problem data to compute j-th
+ solution derivative at k0;
+ - solve: get HF solution as complex numpy vector.
+ HSEngine: Hilbert space general purpose engine. Should have members:
+ - norm: compute Hilbert norm of expression represented as complex
+ numpy vector;
+ - plot: plot Hilbert expression represented as complex numpy vector.
+ k0(optional): Default parameter. Defaults to 0.
+ w(optional): Weight for norm computation. If None, set to Re(k0).
+ Defaults to None.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True;
+ - 'rho': weight for computation of original Pade' approximant;
+ defaults to inf, i.e. fast approximant;
+ - 'M': degree of Pade' approximant numerator; defaults to 0;
+ - 'N': degree of Pade' approximant denominator; defaults to 0;
+ - 'E': total number of derivatives current approximant relies upon;
+ defaults to Emax;
+ - 'Emax': total number of derivatives of solution map to be
+ computed; defaults to E;
+ - 'robustTol': tolerance for robust Pade' denominator management;
+ defaults to 0;
+ - 'sampleType': label of sampling type; available values are:
+ - 'ARNOLDI': orthogonalization of solution derivatives through
+ Arnoldi algorithm;
+ - 'KRYLOV': standard computation of solution derivatives.
+ Defaults to 'KRYLOV'.
+ Defaults to empty dict.
+ plotDer(optional): Whether to plot derivatives of the Helmholtz
+ solution map. Defaults to False.
+ equilibration(optional): Whether to apply equilibration to Gram matrix.
+ Defaults to False.
+ verboseRobust(optional): Verbosity level for robustness-related
+ parameter modifications. Defaults to 1.
+ cleanupParameters(optional): Parameter values for pole cleanup.
+ Available fields are:
+ - 'boolCondition'(bool function handle): if evaluation on pole
+ returns False, pole is removed; defaults to always True;
+ - 'residueCheck'(bool): whether to apply residue check; defaults to
+ False;
+ - 'residueNPoints'(int): number of sample points to be used in
+ residue estimation; defaults to 2;
+ - 'residueRadius'(float): sample radius to be used in residue
+ estimation by Cauchy formula; defaults to 1e-5;
+ - 'residueTol'(float): target tolerance to be used in residue
+ estimation; defaults to 1e-4.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ HSEngine: Hilbert space general purpose engine.
+ solDerivatives: Array whose columns are FE dofs of solution
+ derivatives.
+ k0: Default parameter.
+ w: Weight for norm computation.
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots;
+ - 'rho': weight for computation of original Pade' approximant;
+ defaults to inf, i.e. fast approximant;
+ - 'M': degree of Pade' approximant numerator;
+ - 'N': degree of Pade' approximant denominator;
+ - 'E': total number of derivatives current approximant relies upon;
+ - 'Emax': total number of derivatives of solution map to be
+ computed;
+ - 'robustTol': tolerance for robust Pade' denominator management;
+ - 'sampleType': label of sampling type.
+ M: Numerator degree of approximant.
+ N: Denominator degree of approximant.
+ E: Number of solution derivatives over which current approximant is
+ based upon.
+ Emax: Total number of solution derivatives to be computed.
+ robustTol: tolerance for robust Pade' denominator management.
+ sampleType: Label of sampling type.
+ plotDer: Whether to plot derivatives of the Helmholtz solution map.
+ initialHFData: HF problem initial data.
+ equilibration: Whether to apply equilibration to Gram matrix.
+ verboseRobust: Verbosity level for robustness-related parameter
+ modifications.
+ cleanupParameters; Parameter values for pole cleanup. Available fields
+ are:
+ - 'boolCondition': if evaluation on pole returns False, pole is
+ removed;
+ - 'residueCheck': whether to apply residue check;
+ - 'residueNPoints': number of sample points to be used in residue
+ estimation;
+ - 'residueRadius': sample radius to be used in residue estimation
+ by Cauchy formula;
+ - 'residueTol': target tolerance to be used in residue estimation.
+ energyNormMatrix: Sparse matrix (in CSC format) assoociated to
+ w-weighted H10 norm.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ G: Square Numpy 2D vector of size (N+1) corresponding to Pade'
+ denominator matrix (see paper).
+ equilPowers: Numpy 1D (column) array containing equilibration powers.
+ Does not exist if equilibration is disabled.
+ Q: Numpy 1D vector containing complex coefficients of approximant
+ denominator.
+ P: Numpy 2D vector whose columns are FE dofs of coefficients of
+ approximant numerator.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ """
+
+ extraApproxParameters = ["M", "N", "robustTol", "rho"]
+
+ def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
+ k0 : complex = 0, w : float = None,
+ approxParameters : dict = {}, plotDer : bool = False,
+ equilibration : bool = False, verboseRobust : int = 1,
+ cleanupParameters : dict = {}):
+ self.equilibration = equilibration
+ self.verboseRobust = verboseRobust
+ self.cleanupParameters = cleanupParameters
+ ROMApproximantTaylor.__init__(self, HFEngine, HSEngine, k0, w,
+ approxParameters, plotDer)
+
+ def parameterList(self) -> list:
+ """
+ List containing keys which are allowed in approxParameters.
+
+ Returns:
+ List of strings.
+ """
+ return (ROMApproximantTaylor.parameterList(self)
+ + ROMApproximantTaylorPade.extraApproxParameters)
+
+ @property
+ def approxParameters(self):
+ """
+ Value of approximant parameters. Its assignment may change M, N, E,
+ Emax and robustTol.
+ """
+ return self._approxParameters
+ @approxParameters.setter
+ def approxParameters(self, approxParams):
+ if not hasattr(self, "approxParameters"):
+ self._approxParameters = {}
+ approxParameters = utilities.purgeDict(approxParams,
+ ROMApproximantTaylorPade.parameterList(self),
+ dictname = "approxParameters")
+ keyList = list(approxParameters.keys())
+ if "M" in keyList:
+ self.M = 0 #to avoid warnings if M and E are changed simultaneously
+ if "N" in keyList:
+ self.N = 0 #to avoid warnings if N and E are changed simultaneously
+ approxParametersCopy = utilities.purgeDict(approxParameters,
+ ROMApproximantTaylorPade.extraApproxParameters,
+ True, True)
+ ROMApproximantTaylor.approxParameters.fset(self, approxParametersCopy)
+ if "rho" in keyList:
+ self.rho = approxParameters["rho"]
+ elif hasattr(self, "rho"):
+ self.rho = self.rho
+ else:
+ self.rho = np.inf
+ if "M" in keyList:
+ self.M = approxParameters["M"]
+ elif hasattr(self, "M"):
+ self.M = self.M
+ else:
+ self.M = 0
+ if "N" in keyList:
+ self.N = approxParameters["N"]
+ elif hasattr(self, "N"):
+ self.N = self.N
+ else:
+ self.N = 0
+ if "robustTol" in keyList:
+ self.robustTol = approxParameters["robustTol"]
+ elif hasattr(self, "robustTol"):
+ self.robustTol = self.robustTol
+ else:
+ self.robustTol = 0
+
+ @property
+ def rho(self):
+ """Value of rho."""
+ return self._rho
+ @rho.setter
+ def rho(self, rho):
+ self._rho = np.abs(rho)
+ self._approxParameters["rho"] = self.rho
+
+ @property
+ def M(self):
+ """Value of M. Its assignment may change Emax and E."""
+ return self._M
+ @M.setter
+ def M(self, M):
+ if M < 0: raise ArithmeticError("M must be non-negative.")
+ self._M = M
+ self._approxParameters["M"] = self.M
+ if hasattr(self, "Emax") and self.Emax < self.M:
+ warnings.warn("Prescribed Emax is too small. Updating Emax to M.")
+ self.Emax = self.M
+ if hasattr(self, "E") and self.E < self.M:
+ warnings.warn("Prescribed E is too small. Updating E to M.")
+ self.E = self.M
+
+ @property
+ def N(self):
+ """Value of N. Its assignment may change Emax."""
+ return self._N
+ @N.setter
+ def N(self, N):
+ if N < 0: raise ArithmeticError("N must be non-negative.")
+ self._N = N
+ self._approxParameters["N"] = self.N
+ if hasattr(self, "Emax") and self.Emax < self.N:
+ warnings.warn("Prescribed Emax is too small. Updating Emax to N.")
+ self.Emax = self.N
+ if hasattr(self, "E") and self.E < self.N:
+ warnings.warn("Prescribed E is too small. Updating E to N.")
+ self.E = self.N
+
+ @property
+ def robustTol(self):
+ """Value of tolerance for robust Pade' denominator management."""
+ return self._robustTol
+ @robustTol.setter
+ def robustTol(self, robustTol):
+ if robustTol < 0.:
+ warnings.warn(("Overriding prescribed negative robustness "
+ "tolerance to 0."))
+ robustTol = 0.
+ self._robustTol = robustTol
+ self._approxParameters["robustTol"] = self.robustTol
+
+ @property
+ def E(self):
+ """Value of E. Its assignment may change Emax."""
+ return self._E
+ @E.setter
+ def E(self, E):
+ if E < 0: raise ArithmeticError("E must be non-negative.")
+ self._E = E
+ if hasattr(self, "M") and self.M > self.E:
+ warnings.warn("Prescribed E is too small. Updating E to M.")
+ self.E = self.M
+ if hasattr(self, "N") and self.N > self.E:
+ warnings.warn("Prescribed E is too small. Updating E to N.")
+ self.E = self.N
+ self._approxParameters["E"] = self.E
+ if hasattr(self, "Emax") and self.Emax < self.E:
+ warnings.warn("Prescribed Emax is too small. Updating Emax to E.")
+ self.Emax = self.E
+
+ @property
+ def Emax(self):
+ """Value of Emax. Its assignment may reset computed derivatives."""
+ return self._Emax
+ @Emax.setter
+ def Emax(self, Emax):
+ if Emax < 0: raise ArithmeticError("Emax must be non-negative.")
+ if hasattr(self, "Emax"): EmaxOld = self.Emax
+ else: EmaxOld = -1
+ vals, label = [0] * 3, {0:"E", 1:"M", 2:"N"}
+ if hasattr(self, "E"): vals[0] = self.E
+ if hasattr(self, "M"): vals[1] = self.M
+ if hasattr(self, "N"): vals[2] = self.N
+ idxmax = np.argmax(vals)
+ if vals[idxmax] > Emax:
+ warnings.warn("Prescribed Emax is too small. Updating Emax to {}."\
+ .format(label[idxmax]))
+ self.Emax = vals[idxmax]
+ else:
+ self._Emax = Emax
+ self._approxParameters["Emax"] = Emax
+ if (EmaxOld > self.Emax and self.solDerivatives is not None
+ and self.HArnoldi is not None and self.RArnoldi is not None):
+ self.solDerivatives = self.solDerivatives[:, : self.Emax + 1]
+ if self.sampleType == "ARNOLDI":
+ self.HArnoldi = self.HArnoldi[: self.Emax + 1,
+ : self.Emax + 1]
+ self.RArnoldi = self.RArnoldi[: self.Emax + 1,
+ : self.Emax + 1]
+ else:
+ self.resetSamples()
+
+ @property
+ def cleanupParameters(self):
+ """Value of cleanupParameters."""
+ return self._cleanupParameters
+ @cleanupParameters.setter
+ def cleanupParameters(self, cleanupParameters):
+ allowedCleanupKeys = ['boolCondition','residueCheck','residueNPoints',
+ 'residueRadius','residueTol']
+ cleanupKeys = cleanupParameters.keys()
+ if 'boolCondition' not in cleanupKeys:
+ cleanupParameters['boolCondition'] = lambda x : True
+ if 'residueCheck' not in cleanupKeys:
+ cleanupParameters['residueCheck'] = False
+ if 'residueNPoints' not in cleanupKeys:
+ cleanupParameters['residueNPoints'] = 2
+ if 'residueRadius' not in cleanupKeys:
+ cleanupParameters['residueRadius'] = 1e-5
+ if 'residueTol' not in cleanupKeys:
+ cleanupParameters['residueTol'] = 1e-4
+ cleanupParameters['boolCondition'] = np.vectorize(
+ cleanupParameters['boolCondition'])
+ self._cleanupParameters = {key : cleanupParameters[key]
+ for key in allowedCleanupKeys}
+
+ def setupApprox(self):
+ """
+ Compute Pade' approximant.
+ See Section 4 in Fast Pade' paper.
+ SVD-based robust eigenvalue management.
+ """
+ self.computeDerivatives()
+ if not self.checkComputedApprox():
+ while True:
+ if self.POD:
+ if self.sampleType == "ARNOLDI":
+ ev, eV = self.findeveVGArnoldi()
+ else:
+ ev, eV = self.findeveVGQR()
+ else:
+ self.buildG()
+ ev, eV = np.linalg.eigh(self.G)
+ ts = self.robustTol * np.linalg.norm(ev)
+ Nnew = np.sum(np.abs(ev) >= ts)
+ diff = self.N - Nnew
+ if diff <= 0: break
+ Enew = self.E - diff
+ Mnew = min(self.M, Enew)
+ if Mnew == self.M:
+ strM = ""
+ else:
+ strM = ", M from {0} to {1},".format(self.M, Mnew)
+ if self.verboseRobust >= 1:
+ print(("Smallest {0} eigenvalues below tolerance.\n"
+ "Reducing N from {1} to {2}{5} and E from {3} to "
+ "{4}.").format(diff + 1, self.N, Nnew,
+ self.E, Enew, strM))
+ newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew}
+ self.approxParameters = newParameters
+
+
+ self.Q = eV[::-1, 0]
+ if self.equilibration:
+ self.Q = np.multiply(self.equilPowers.flatten()[::-1], self.Q)
+
+ self._cleanup()
+
+ QToeplitz = np.zeros((self.Emax + 1, self.M + 1),
+ dtype = np.complex)
+ for i in range(len(self.Q)):
+ rng = np.arange(self.M + 1 - i)
+ QToeplitz[rng, rng + i] = self.Q[i]
+ if self.sampleType == "ARNOLDI":
+ QToeplitz = self.RArnoldi.dot(QToeplitz)
+ self.P = self.solDerivatives.dot(QToeplitz)
+
+ self.lastApproxParameters = copy(self.approxParameters)
+
+ def buildG(self):
+ """Assemble Pade' denominator matrix."""
+ if self.N == 0:
+ self.G = np.array([[1]], dtype = np.complex)
+ return
+ self.computeDerivatives()
+ if np.isinf(self.rho):
+ lgSize = self.N + 1
+ else:
+ lgSize = self.E + self.N - self.M
+ Gshift = self.E - lgSize + 1
+
+ largeG = np.zeros((lgSize, lgSize), dtype=np.complex)
+ largeG = self.solDerivatives[:, Gshift:self.E+1].conj().T.dot(
+ self.energyNormMatrix.dot(self.solDerivatives[:, Gshift:self.E+1]))
+ if np.isinf(self.rho):
+ self.G = largeG
+ else:
+ self.G = np.zeros((self.N + 1, self.N + 1), dtype=np.complex)
+ for k in range(self.E - self.M):
+ self.G = (self.G + self.rho ** (2 * k)
+ * largeG[k : k+self.N+1, k : k+self.N+1])
+ if self.equilibration:
+ Gd = np.diag(self.G)
+ gamma = np.average(np.abs(np.divide(Gd[1:], Gd[:-1])),
+ weights = np.power(1.2, np.arange(self.N)))**.5
+ self.equilPowers = np.power(gamma, np.arange(self.N + 1)[:, None])
+ self.G = np.multiply(self.equilPowers,
+ np.multiply(self.equilPowers, self.G).T).T
+
+ def _cleanup(self):
+ """Cleanup Pade' denominator by removing unwanted poles."""
+ if self.N == 0: return
+ poles = np.roots(self.Q[::-1]) + self.k0
+ NpolesOld = len(poles)
+ poles = poles[self.cleanupParameters['boolCondition'](poles)]
+
+ if self.cleanupParameters['residueCheck']:
+ resR = self.cleanupParameters['residueRadius']
+ resTol = self.cleanupParameters['residueTol']
+ residues = np.zeros_like(poles)
+ for l, pole in enumerate(poles):
+ resV = np.zeros(self.solDerivatives.shape[0],
+ dtype = np.complex)
+ NPoints = self.cleanupParameters['residueNPoints']
+ for theta in 2 * PI * np.arange(NPoints) / NPoints:
+ deltapole = resR * np.exp(1.j * theta)
+ ksample = pole + deltapole
+ self.solveHF(ksample)
+ resV = deltapole / NPoints * (resV + self.uHF)
+ residues[l] = np.abs(resV.dot(
+ self.energyNormMatrix.dot(resV).conj())) ** .5
+ poles = poles[residues >= resTol]
+
+ NpolesNew = len(poles)
+ if NpolesOld > NpolesNew:
+ if self.verboseRobust >= 1:
+ sSing = "s" * (NpolesOld - NpolesNew > 1)
+ print(("Identified {0} pole{1} to be removed.\n"
+ "Reducing N from {2} to {3}.").format(
+ NpolesOld - NpolesNew, sSing, NpolesOld, NpolesNew))
+ self.N = NpolesNew
+ newQ = np.polyfit(np.append(poles - self.k0, 0.),
+ np.append(np.zeros(NpolesNew), 1.), NpolesNew)
+ self.Q = newQ[::-1] / np.linalg.norm(newQ)
+
+ def findeveVGQR(self):
+ """
+ Compute eigenvalues and eigenvectors of Pade' denominator matrix
+ through SVD of R factor. See ``Householder triangularization of a
+ quasimatrix'', L.Trefethen, 2008 for QR algorithm.
+
+ Returns:
+ Eigenvalues in ascending order and corresponding eigenvector
+ matrix.
+ """
+ self.computeDerivatives()
+ A = copy(self.solDerivatives[:, (self.E - self.N) : (self.E + 1)])
+ M = self.energyNormMatrix
+ E = np.zeros_like(A, dtype = np.complex)
+ R = np.zeros((self.N + 1, self.N + 1), dtype = np.complex)
+ for k in range(self.N + 1):
+ E[:, k] = (np.random.randn(E.shape[0])
+ + 1.j * np.random.randn(E.shape[0]))
+ if k > 0:
+ for times in range(2):
+ E[:, k] = E[:, k] - E[:, :k].dot(
+ (E[:, k].conj().dot(M.dot(E[:, :k]))).conj())
+ E[:, k] = E[:, k] / (E[:, k].conj().dot(M.dot(E[:, k]))) ** .5
+ R[k, k] = np.abs(A[:, k].conj().dot(M.dot(A[:, k]))) ** .5
+ alpha = E[:, k].conj().dot(M.dot(A[:, k]))
+ if np.isclose(np.abs(alpha), 0.): s = 1.
+ else: s = - alpha / np.abs(alpha)
+ E[:, k] = s * E[:, k]
+ v = R[k, k] * E[:, k] - A[:, k]
+ for times in range(2):
+ v = v - E[:, :k].dot((v.conj().dot(M.dot(E[:, :k]))).conj())
+ sigma = np.abs(v.conj().dot(M.dot(v))) ** .5
+ if np.isclose(sigma, 0.): v = E[:, k]
+ else: v = v / sigma
+ J = np.arange(k + 1, self.N + 1)
+ vtQ = v.conj().dot(M.dot(A[:, J]))
+ A[:, J] = A[:, J] - 2 * np.outer(v, vtQ)
+ R[k, J] = E[:, k].conj().dot(M.dot(A[:, J]))
+ A[:, J] = A[:, J] - np.outer(E[:, k], R[k, J])
+ if self.equilibration:
+ Rd = np.diag(R)
+ gamma = np.average(np.abs(np.divide(Rd[1:], Rd[:-1])),
+ weights = np.power(1.2, np.arange(self.N)))**.5
+ self.equilPowers = np.power(gamma, np.arange(self.N + 1)[:, None])
+ R = np.multiply(self.equilPowers, R.T).T
+ _, s, V = np.linalg.svd(R, full_matrices = False)
+ return s[::-1], V.conj().T[:, ::-1]
+
+ def findeveVGArnoldi(self):
+ """
+ Compute eigenvalues and eigenvectors of Pade' denominator matrix
+ through SVD of R factor of solution derivatives orthogonalized by
+ Arnoldi algorithm.
+
+ Returns:
+ Eigenvalues in ascending order and corresponding eigenvector
+ matrix.
+ """
+ self.computeDerivatives()
+ if self.sampleType != "ARNOLDI":
+ raise Exception(("Eigensolver different from 'ARNOLDI'."
+ "Arnoldi eigenproblem solver cannot be called."))
+ R = self.RArnoldi[: self.E + 1, self.E - self.N : self.E + 1]
+ if self.equilibration:
+ Rd = np.diag(R)
+ gamma = np.average(np.abs(np.divide(Rd[1:], Rd[:-1])),
+ weights = np.power(1.2, np.arange(self.N)))**.5
+ self.equilPowers = np.power(gamma, np.arange(self.N + 1)[:, None])
+ R = np.multiply(self.equilPowers, R.T).T
+ _, s, V = np.linalg.svd(R, full_matrices = False)
+ return s[::-1], V.conj().T[:, ::-1]
+
+ def evalApprox(self, k:complex) -> ("Fenics function", "Fenics function"):
+ """
+ Evaluate Pade' approximant at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ Approximant as numpy complex vector.
+ """
+ self.setupApprox()
+ powerlist = np.power(k - self.k0, range(max(self.M, self.N) + 1))
+ self.uApp = (self.P.dot(powerlist[:self.M + 1])
+ / self.Q.dot(powerlist[:self.N + 1]))
+ return self.uApp
+
+ def getPoles(self, centered : bool = False) -> "numpy 1D array":
+ """
+ Obtain approximant poles.
+
+ Args:
+ centered(optional): Whether to return pole positions relative to
+ approximation center. Defaults to False.
+
+ Returns:
+ Numpy complex vector of poles.
+ """
+ self.setupApprox()
+ return np.roots(self.Q[::-1]) + self.k0 * centered
+
+
\ No newline at end of file
diff --git a/main/ROMApproximantTaylorRB.py b/main/ROMApproximantTaylorRB.py
new file mode 100644
index 0000000..9a48054
--- /dev/null
+++ b/main/ROMApproximantTaylorRB.py
@@ -0,0 +1,303 @@
+#!/usr/bin/python
+
+from copy import copy
+import warnings
+import numpy as np
+import scipy as sp
+import utilities
+from ROMApproximantTaylor import ROMApproximantTaylor
+
+class ROMApproximantTaylorRB(ROMApproximantTaylor):
+ """
+ ROM single-point fast RB approximant computation for parametric problems
+ with polynomial dependence up to degree 2.
+
+ Args:
+ HFEngine: HF problem solver. Should have members:
+ - energyNormMatrix: sparse matrix (in CSC format) associated to
+ w-weighted H10 norm;
+ - problemData: list of HF problem data (excluding k);
+ - setProblemData: set HF problem data and k to prescribed values;
+ - getLSBlocks: get blocks of HF linear system as sparse matrices in
+ CSC format;
+ - liftDirichletData: perform lifting of Dirichlet BC as numpy
+ complex vector;
+ - setupDerivativeComputation: setup HF problem data to compute j-th
+ solution derivative at k0;
+ - solve: get HF solution as complex numpy vector.
+ HSEngine: Hilbert space general purpose engine. Should have members:
+ - norm: compute Hilbert norm of expression represented as complex
+ numpy vector;
+ - plot: plot Hilbert expression represented as complex numpy vector.
+ k0(optional): Default parameter. Defaults to 0.
+ w(optional): Weight for norm computation. If None, set to Re(k0).
+ Defaults to None.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True;
+ - 'R': rank for Galerkin projection; defaults to E + 1;
+ - 'E': total number of derivatives current approximant relies upon;
+ defaults to Emax;
+ - 'Emax': total number of derivatives of solution map to be
+ computed; defaults to E;
+ - 'sampleType': label of sampling type; available values are:
+ - 'ARNOLDI': orthogonalization of solution derivatives through
+ Arnoldi algorithm;
+ - 'KRYLOV': standard computation of solution derivatives.
+ Defaults to 'KRYLOV'.
+ Defaults to empty dict.
+ plotDer(optional): Whether to plot derivatives of the Helmholtz
+ solution map. Defaults to False.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ HSEngine: Hilbert space general purpose engine.
+ solDerivatives: Array whose columns are FE dofs of solution
+ derivatives.
+ k0: Default parameter.
+ w: Weight for norm computation.
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots;
+ - 'R': rank for Galerkin projection;
+ - 'E': total number of derivatives current approximant relies upon;
+ - 'Emax': total number of derivatives of solution map to be
+ computed;
+ - 'sampleType': label of sampling type.
+ R: Rank for Galerkin projection.
+ E: Number of solution derivatives over which current approximant is
+ based upon.
+ Emax: Total number of solution derivatives to be computed.
+ sampleType: Label of sampling type, i.e. 'KRYLOV'.
+ plotDer: Whether to plot derivatives of the Helmholtz solution map.
+ energyNormMatrix: Sparse matrix (in CSC format) associated to
+ w-weighted H10 norm.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ solLifting: Numpy complex vector with lifting of real part of Dirichlet
+ boundary datum.
+ projMat: Numpy matrix representing projection onto RB space.
+ projMat: Numpy matrix representing projection onto RB space.
+ As: List of sparse matrices (in CSC format) representing coefficients
+ of linear system matrix wrt k.
+ bs: List of numpy vectors representing coefficients of linear system
+ RHS wrt k.
+ ARBs: List of sparse matrices (in CSC format) representing RB
+ coefficients of linear system matrix wrt k.
+ bRBs: List of numpy vectors representing RB coefficients of linear
+ system RHS wrt k.
+ """
+
+ extraApproxParameters = ["R"]
+
+ def __init__(self, HFEngine:'HF solver', HSEngine:'HS engine',
+ k0 : complex = 0, w : float = None,
+ approxParameters : dict = {}, plotDer : bool = False):
+ ROMApproximantTaylor.__init__(self, HFEngine, HSEngine, k0, w,
+ approxParameters, plotDer)
+ self.solLifting = self.HFEngine.liftDirichletData()
+
+ def resetSamples(self):
+ """Reset samples."""
+ ROMApproximantTaylor.resetSamples(self)
+ self.projMat = None
+
+ def parameterList(self) -> list:
+ """
+ List containing keys which are allowed in approxParameters.
+
+ Returns:
+ List of strings.
+ """
+ return (ROMApproximantTaylor.parameterList(self)
+ + ROMApproximantTaylorRB.extraApproxParameters)
+
+ @property
+ def approxParameters(self):
+ """
+ Value of approximant parameters. Its assignment may change M, N and S.
+ """
+ return self._approxParameters
+ @approxParameters.setter
+ def approxParameters(self, approxParams):
+ if not hasattr(self, "approxParameters"):
+ self._approxParameters = {}
+ approxParameters = utilities.purgeDict(approxParams,
+ ROMApproximantTaylorRB.parameterList(self),
+ dictname = "approxParameters")
+ keyList = list(approxParameters.keys())
+ approxParametersCopy = utilities.purgeDict(approxParameters,
+ ROMApproximantTaylorRB.extraApproxParameters,
+ True, True)
+ ROMApproximantTaylor.approxParameters.fset(self, approxParametersCopy)
+ if "R" in keyList:
+ self.R = approxParameters["R"]
+ elif hasattr(self, "R"):
+ self.R = self.R
+ else:
+ self.R = self.E + 1
+
+ @property
+ def R(self):
+ """Value of R. Its assignment may change S."""
+ return self._R
+ @R.setter
+ def R(self, R):
+ if R < 0: raise ArithmeticError("R must be non-negative.")
+ self._R = R
+ self._approxParameters["R"] = self.R
+ if hasattr(self, "E") and self.E + 1 < self.R:
+ warnings.warn("Prescribed E is too small. Updating E to R - 1.")
+ self.E = self.R - 1
+
+ def manageDerivatives(self, u:"2-tuple of Fenics function", pos:int)\
+ -> ("Fenics function", "Fenics function"):
+ """
+ Store derivatives of solution map. Subtract lifting for order 0.
+
+ Args:
+ u: 2-tuple containing real and imaginary parts of FE dofs of
+ derivative.
+ pos: Derivative index.
+
+ Returns:
+ Real and imaginary parts of derivative (possibly adjusted).
+ """
+ if pos == 0:
+ self.As, self.bs = self.HFEngine.getLSBlocks()
+ u = u - self.solLifting
+ return ROMApproximantTaylor.manageDerivatives(self, u, pos)
+
+ def setupApprox(self):
+ """
+ Setup RB system. For usage of correlation matrix in POD see Section
+ 6.3.1 in 'Reduced Basis Methods for PDEs. An introduction', A.
+ Quarteroni, A. Manzoni, F. Negri, Springer, 2016.
+ """
+ need2Setup = (self.solDerivatives is None) or (self.projMat is None)
+ self.computeDerivatives()
+ if need2Setup:
+ if self.POD and not self.sampleType == "ARNOLDI":
+ A = copy(self.solDerivatives)
+ M = self.energyNormMatrix
+ V = np.zeros_like(A, dtype = np.complex)
+ Q = np.zeros_like(A, dtype = np.complex)
+ R = np.zeros((A.shape[1], A.shape[1]), dtype = np.complex)
+ for k in range(A.shape[1]):
+ Q[:, k] = (np.random.randn(Q.shape[0])
+ + 1.j * np.random.randn(Q.shape[0]))
+ if k > 0:
+ for times in range(2):
+ Q[:, k] = Q[:, k] - Q[:, :k].dot(
+ (Q[:, k].conj().dot(M.dot(Q[:, :k]))).conj())
+ Q[:, k] = Q[:, k]/(Q[:, k].conj().dot(M.dot(Q[:, k])))**.5
+ R[k, k] = np.abs(A[:, k].conj().dot(M.dot(A[:, k]))) ** .5
+ alpha = Q[:, k].conj().dot(M.dot(A[:, k]))
+ if np.isclose(np.abs(alpha), 0.): s = 1.
+ else: s = - alpha / np.abs(alpha)
+ Q[:, k] = s * Q[:, k]
+ v = R[k, k] * Q[:, k] - A[:, k]
+ for times in range(2):
+ v = v - Q[:, :k].dot((v.conj().dot(M.dot(Q[:, :k]))
+ ).conj())
+ sigma = np.abs(v.conj().dot(M.dot(v))) ** .5
+ if np.isclose(sigma, 0.): v = Q[:, k]
+ else: v = v / sigma
+ V[:, k] = v
+ J = np.arange(k + 1, A.shape[1])
+ vtQ = v.conj().dot(M.dot(A[:, J]))
+ A[:, J] = A[:, J] - 2 * np.outer(v, vtQ)
+ R[k, J] = Q[:, k].conj().dot(M.dot(A[:, J]))
+ A[:, J] = A[:, J] - np.outer(Q[:, k], R[k, J])
+ for k in range(A.shape[1] - 1, -1, -1):
+ v = V[:, k]
+ J = np.arange(k, A.shape[1])
+ vtQ = v.conj().dot(M.dot(Q[:, J]))
+ Q[:, J] = Q[:, J] - 2 * np.outer(v, vtQ)
+ self.projMatQ = Q
+ self.projMatR = R
+ if self.POD and not self.sampleType == "ARNOLDI":
+ U, _, _ = np.linalg.svd(self.projMatR[: self.R, : self.R])
+ self.projMat = self.projMatQ[:, : self.R].dot(U)
+ else:
+ self.projMat = self.solDerivatives[:, : self.R]
+ self.assembleReducedSystem()
+ self.lastApproxParameters = copy(self.approxParameters)
+
+ def assembleReducedSystem(self):
+ """Build affine blocks of RB linear system through projections."""
+ projMatH = self.projMat.T.conjugate()
+ self.ARBs = [None] * len(self.As)
+ self.bRBs = [None] * max(len(self.As), len(self.bs))
+ for j in range(len(self.As)):
+ self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat))
+ if j < len(self.bs):
+ self.bRBs[j] = projMatH.dot(self.bs[j]
+ - self.As[j].dot(self.solLifting))
+ else:
+ self.bRBs[j] = - projMatH.dot(self.As[j].dot(self.solLifting))
+ for j in range(len(self.As), len(self.bs)):
+ self.bRBs[j] = projMatH.dot(self.bs[j])
+
+ def solveReducedSystem(self, k:complex) -> "Numpy 1D array":
+ """
+ Solve RB linear system.
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ Solution of RB linear system.
+ """
+ self.setupApprox()
+ ARBk = self.ARBs[0][: self.R, : self.R]
+ bRBk = self.bRBs[0][: self.R]
+ for j in range(1, len(self.ARBs)):
+ ARBk = ARBk + np.power(k, j) * self.ARBs[j][:self.R, :self.R]
+ for j in range(1, len(self.bRBs)):
+ bRBk = bRBk + np.power(k, j) * self.bRBs[j][:self.R]
+ return self.projMat[:, :self.R].dot(np.linalg.solve(ARBk, bRBk))
+
+ def evalApprox(self, k:complex) -> ("Fenics function", "Fenics function"):
+ """
+ Evaluate RB approximant at arbitrary wavenumber.
+
+ Args:
+ k: Target wavenumber.
+
+ Returns:
+ Real and imaginary parts of approximant.
+ """
+ self.setupApprox()
+ self.uApp = self.solLifting + self.solveReducedSystem(k)
+ return self.uApp
+
+ def getPoles(self, centered : bool = False) -> "numpy 1D array":
+ """
+ Obtain approximant poles.
+
+ Args:
+ centered(optional): Whether to return pole positions relative to
+ approximation center. Defaults to False.
+
+ Returns:
+ Numpy complex vector of poles.
+ """
+ self.setupApprox()
+ A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1),
+ dtype = np.complex)
+ B = np.zeros_like(A)
+ A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0]
+ for j in range(len(self.ARBs) - 1):
+ Aj = self.ARBs[j + 1]
+ B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj
+ II = np.arange(self.ARBs[0].shape[0],
+ self.ARBs[0].shape[0] * (len(self.ARBs) - 1))
+ B[II, II - self.ARBs[0].shape[0]] = 1.
+ return sp.linalg.eigvals(A, B) - self.k0 * (not centered)
+
diff --git a/main/__init__.py b/main/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/main/utilities.py b/main/utilities.py
new file mode 100644
index 0000000..f405d07
--- /dev/null
+++ b/main/utilities.py
@@ -0,0 +1,90 @@
+#!/usr/bin/python
+
+import warnings
+import numpy as np
+#from copy import copy
+
+def findDictStrKey(key, keyList):
+ for akey in keyList:
+ if isinstance(key, str) and key.lower() == akey.lower():
+ return akey
+ return None
+
+def purgeList(lst:list, allowedEntries : list = [], silent : bool = False,
+ complement : bool = False, listname : str = ""):
+ if listname != "":
+ listname = " in " + listname
+ allowedEntriesSet = frozenset(allowedEntries)
+ lstcp = []
+ for x in lst:
+ if (x in allowedEntriesSet) != complement:
+ lstcp = lstcp + [x]
+ elif not silent:
+ warnings.warn("Ignoring entry {0}{1}.".format(x, listname))
+ return lstcp
+
+def purgeDict(dct:dict, allowedKeys : list = [], silent : bool = False,
+ complement : bool = False, dictname : str = ""):
+ if dictname != "":
+ dictname = " in " + dictname
+ dctcp = {}
+ for key in dct.keys():
+ akey = findDictStrKey(key, allowedKeys)
+ if (akey is None) != complement:
+ if not silent:
+ warnings.warn("Ignoring key {0}{2} with value {1}."\
+ .format(key, dct[key], dictname))
+ else:
+ if akey is None:
+ akey = key
+ dctcp[akey] = dct[key]
+ return dctcp
+
+prime_v = [] #memoization vector
+
+def squareResonances(a:int, b:int, zero_terms : bool = True):
+ spectrum = []
+ a = max(int(np.floor(a)), 0)
+ b = max(int(np.ceil(b)), 0)
+ global prime_v
+ if len(prime_v) == 0:
+ prime_v = [2, 3]
+ if a > prime_v[-1]:
+ for i in range(prime_v[-1], a, 2):
+ get_next_prime_factor(i)
+ for i in range(a, b + 1):
+ spectrum = spectrum + [i] * count_square_sums(i, zero_terms)
+ return np.array(spectrum)
+
+def get_next_prime_factor(n):
+ global prime_v
+ for x in prime_v:
+ if n % x == 0:
+ return x
+ if prime_v[-1] < n:
+ prime_v = prime_v + [n]
+ return n
+
+def prime_factorize(n):
+ factors = []
+ number = n
+ while number > 1:
+ factor = get_next_prime_factor(number)
+ factors.append(factor)
+ number = int(number / factor)
+ if n < -1:
+ factors[0] = -factors[0]
+ return list(factors)
+
+def count_square_sums(n, zero_terms : bool = True):
+ if n < 2: return (n + 1) * zero_terms
+ factors = prime_factorize(n)
+ funique, fcounts = np.unique(factors, return_counts = True)
+ count = 1
+ for fac, con in zip(funique, fcounts): #using number theory magic
+ if fac % 4 == 1:
+ count = count * (con + 1)
+ elif fac % 4 == 3 and con % 2 == 1:
+ return 0
+ return count + (2 * zero_terms - 1) * (int(n ** .5) ** 2 == n)
+