pax_global_header00006660000000000000000000000064131562277650014530gustar00rootroot0000000000000052 comment=626894c5bf93127a94db1ce737abd754cc59903f
cim-master/000077500000000000000000000000001315622776500131355ustar00rootroot00000000000000cim-master/LICENSE.txt000066400000000000000000000432541315622776500147700ustar00rootroot00000000000000 GNU GENERAL PUBLIC LICENSE
Version 2, June 1991
Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
Preamble
The licenses for most software are designed to take away your
freedom to share and change it. By contrast, the GNU General Public
License is intended to guarantee your freedom to share and change free
software--to make sure the software is free for all its users. This
General Public License applies to most of the Free Software
Foundation's software and to any other program whose authors commit to
using it. (Some other Free Software Foundation software is covered by
the GNU Lesser General Public License instead.) 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
this service 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 make restrictions that forbid
anyone to deny you these rights or to ask you to surrender the rights.
These restrictions translate to certain responsibilities for you if you
distribute copies of the software, or if you modify it.
For example, if you distribute copies of such a program, whether
gratis or for a fee, you must give the recipients all the rights that
you have. 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.
We protect your rights with two steps: (1) copyright the software, and
(2) offer you this license which gives you legal permission to copy,
distribute and/or modify the software.
Also, for each author's protection and ours, we want to make certain
that everyone understands that there is no warranty for this free
software. If the software is modified by someone else and passed on, we
want its recipients to know that what they have is not the original, so
that any problems introduced by others will not reflect on the original
authors' reputations.
Finally, any free program is threatened constantly by software
patents. We wish to avoid the danger that redistributors of a free
program will individually obtain patent licenses, in effect making the
program proprietary. To prevent this, we have made it clear that any
patent must be licensed for everyone's free use or not licensed at all.
The precise terms and conditions for copying, distribution and
modification follow.
GNU GENERAL PUBLIC LICENSE
TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
0. This License applies to any program or other work which contains
a notice placed by the copyright holder saying it may be distributed
under the terms of this General Public License. The "Program", below,
refers to any such program or work, and a "work based on the Program"
means either the Program or any derivative work under copyright law:
that is to say, a work containing the Program or a portion of it,
either verbatim or with modifications and/or translated into another
language. (Hereinafter, translation is included without limitation in
the term "modification".) Each licensee is addressed as "you".
Activities other than copying, distribution and modification are not
covered by this License; they are outside its scope. The act of
running the Program is not restricted, and the output from the Program
is covered only if its contents constitute a work based on the
Program (independent of having been made by running the Program).
Whether that is true depends on what the Program does.
1. You may copy and distribute 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 and disclaimer of warranty; keep intact all the
notices that refer to this License and to the absence of any warranty;
and give any other recipients of the Program a copy of this License
along with the Program.
You may charge a fee for the physical act of transferring a copy, and
you may at your option offer warranty protection in exchange for a fee.
2. You may modify your copy or copies of the Program or any portion
of it, thus forming a work based on the Program, and copy and
distribute such modifications or work under the terms of Section 1
above, provided that you also meet all of these conditions:
a) You must cause the modified files to carry prominent notices
stating that you changed the files and the date of any change.
b) You must cause any work that you distribute or publish, that in
whole or in part contains or is derived from the Program or any
part thereof, to be licensed as a whole at no charge to all third
parties under the terms of this License.
c) If the modified program normally reads commands interactively
when run, you must cause it, when started running for such
interactive use in the most ordinary way, to print or display an
announcement including an appropriate copyright notice and a
notice that there is no warranty (or else, saying that you provide
a warranty) and that users may redistribute the program under
these conditions, and telling the user how to view a copy of this
License. (Exception: if the Program itself is interactive but
does not normally print such an announcement, your work based on
the Program is not required to print an announcement.)
These requirements apply to the modified work as a whole. If
identifiable sections of that work are not derived from the Program,
and can be reasonably considered independent and separate works in
themselves, then this License, and its terms, do not apply to those
sections when you distribute them as separate works. But when you
distribute the same sections as part of a whole which is a work based
on the Program, the distribution of the whole must be on the terms of
this License, whose permissions for other licensees extend to the
entire whole, and thus to each and every part regardless of who wrote it.
Thus, it is not the intent of this section to claim rights or contest
your rights to work written entirely by you; rather, the intent is to
exercise the right to control the distribution of derivative or
collective works based on the Program.
In addition, mere aggregation of another work not based on the Program
with the Program (or with a work based on the Program) on a volume of
a storage or distribution medium does not bring the other work under
the scope of this License.
3. You may copy and distribute the Program (or a work based on it,
under Section 2) in object code or executable form under the terms of
Sections 1 and 2 above provided that you also do one of the following:
a) Accompany it with the complete corresponding machine-readable
source code, which must be distributed under the terms of Sections
1 and 2 above on a medium customarily used for software interchange; or,
b) Accompany it with a written offer, valid for at least three
years, to give any third party, for a charge no more than your
cost of physically performing source distribution, a complete
machine-readable copy of the corresponding source code, to be
distributed under the terms of Sections 1 and 2 above on a medium
customarily used for software interchange; or,
c) Accompany it with the information you received as to the offer
to distribute corresponding source code. (This alternative is
allowed only for noncommercial distribution and only if you
received the program in object code or executable form with such
an offer, in accord with Subsection b above.)
The source code for a work means the preferred form of the work for
making modifications to it. For an executable work, complete source
code means all the source code for all modules it contains, plus any
associated interface definition files, plus the scripts used to
control compilation and installation of the executable. However, as a
special exception, the source code distributed need not include
anything that is normally distributed (in either source or binary
form) with the major components (compiler, kernel, and so on) of the
operating system on which the executable runs, unless that component
itself accompanies the executable.
If distribution of executable or object code is made by offering
access to copy from a designated place, then offering equivalent
access to copy the source code from the same place counts as
distribution of the source code, even though third parties are not
compelled to copy the source along with the object code.
4. You may not copy, modify, sublicense, or distribute the Program
except as expressly provided under this License. Any attempt
otherwise to copy, modify, sublicense or distribute the Program is
void, and will automatically terminate your rights under this License.
However, parties who have received copies, or rights, from you under
this License will not have their licenses terminated so long as such
parties remain in full compliance.
5. You are not required to accept this License, since you have not
signed it. However, nothing else grants you permission to modify or
distribute the Program or its derivative works. These actions are
prohibited by law if you do not accept this License. Therefore, by
modifying or distributing the Program (or any work based on the
Program), you indicate your acceptance of this License to do so, and
all its terms and conditions for copying, distributing or modifying
the Program or works based on it.
6. Each time you redistribute the Program (or any work based on the
Program), the recipient automatically receives a license from the
original licensor to copy, distribute or modify the Program subject to
these terms and conditions. You may not impose any further
restrictions on the recipients' exercise of the rights granted herein.
You are not responsible for enforcing compliance by third parties to
this License.
7. If, as a consequence of a court judgment or allegation of patent
infringement or for any other reason (not limited to patent issues),
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
distribute so as to satisfy simultaneously your obligations under this
License and any other pertinent obligations, then as a consequence you
may not distribute the Program at all. For example, if a patent
license would not permit royalty-free redistribution of the Program by
all those who receive copies directly or indirectly through you, then
the only way you could satisfy both it and this License would be to
refrain entirely from distribution of the Program.
If any portion of this section is held invalid or unenforceable under
any particular circumstance, the balance of the section is intended to
apply and the section as a whole is intended to apply in other
circumstances.
It is not the purpose of this section to induce you to infringe any
patents or other property right claims or to contest validity of any
such claims; this section has the sole purpose of protecting the
integrity of the free software distribution system, which is
implemented by public license practices. Many people have made
generous contributions to the wide range of software distributed
through that system in reliance on consistent application of that
system; it is up to the author/donor to decide if he or she is willing
to distribute software through any other system and a licensee cannot
impose that choice.
This section is intended to make thoroughly clear what is believed to
be a consequence of the rest of this License.
8. If the distribution and/or use of the Program is restricted in
certain countries either by patents or by copyrighted interfaces, the
original copyright holder who places the Program under this License
may add an explicit geographical distribution limitation excluding
those countries, so that distribution is permitted only in or among
countries not thus excluded. In such case, this License incorporates
the limitation as if written in the body of this License.
9. The Free Software Foundation may publish revised and/or new versions
of the 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 a version number of this License which applies to it and "any
later version", you have the option of following the terms and conditions
either of that version or of any later version published by the Free
Software Foundation. If the Program does not specify a version number of
this License, you may choose any version ever published by the Free Software
Foundation.
10. If you wish to incorporate parts of the Program into other free
programs whose distribution conditions are different, write to the author
to ask for permission. For software which is copyrighted by the Free
Software Foundation, write to the Free Software Foundation; we sometimes
make exceptions for this. Our decision will be guided by the two goals
of preserving the free status of all derivatives of our free software and
of promoting the sharing and reuse of software generally.
NO WARRANTY
11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, 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.
12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
REDISTRIBUTE 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.
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
convey 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 2 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, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
Also add information on how to contact you by electronic and paper mail.
If the program is interactive, make it output a short notice like this
when it starts in an interactive mode:
Gnomovision version 69, Copyright (C) year name of author
Gnomovision 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, the commands you use may
be called something other than `show w' and `show c'; they could even be
mouse-clicks or menu items--whatever suits your program.
You should also get your employer (if you work as a programmer) or your
school, if any, to sign a "copyright disclaimer" for the program, if
necessary. Here is a sample; alter the names:
Yoyodyne, Inc., hereby disclaims all copyright interest in the program
`Gnomovision' (which makes passes at compilers) written by James Hacker.
, 1 April 1989
Ty Coon, President of Vice
This 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.
cim-master/README.md000066400000000000000000000107001315622776500144120ustar00rootroot00000000000000cim.py
======
cim.py is a non-linear eigenvalue solver using the contour integral method
proposed by W.-J. Beyn coupled with a [GetDP](http://getdp.info) formulation.
Usage
-----
cim.py pro mesh resolution post origin radius
with:
pro GetDP .pro file
mesh Gmsh .msh file
resolution resolution from .pro file
origin circular contour origin (complex)
radius circular contour radius
Additional optional argument are available. Please run:
cim.py -h
for more information.
Driving [GetDP](http://getdp.info) from cim.py
----------------------------------------------
cim.py and [GetDP](http://getdp.info) are exchanging the following parameters:
angularFreqRe the real part of the frequency
angularFreqIm the imaginary part of the frequency
nRHS the number of right hand side that should be considered (default: 1)
b_i()/b~{i}() the ith right hand side (RHS) of the linear system (first index is 0)
x_i()/x~{i}() the ith solution vector of the linear system (first index is 0)
doInit flag signalling if an initialisation should be done (default: 0)
doSolve flag signalling if a solve should be done for all RHS (default: 0)
doPostpro flag signalling if a post-processing should be done (default: 0)
doApply flag signalling if x() should be applied to the linear system (default: 0)
fileName name of the post-processing file
More information can be found in [bin/solver.py](bin/solver.py). Furthermore,
examples are provided in [example/maxwell\_sibc](example/maxwell_sibc) and
[example/maxwell\_linear](example/maxwell\_linear).
Dependencies
------------
cim.py is a [Python 2.7](https://www.python.org/download/releases/2.7) script.
It relies on the following modules:
argparse
numpy
getdp
The getdp module can be obtained by compiling [GetDP](http://getdp.info).
Among other (and possibly optional) dependencies, [GetDP](http://getdp.info)
relies on [PETSc](https://www.mcs.anl.gov/petsc).
The shell script [dependencies.sh](dependencies.sh) should handle this
automatically if you run:
./dependencies.sh
You will need `wget`, `tar`, `gcc`, `g++` and `gfortran` for this script
to work. Don't forget to update your `PATH`, `PYTHONPATH` and `LD_LIBRARY_PATH`
as instructed.
If you want to handle the [PETSc](https://www.mcs.anl.gov/petsc) and [GetDP](
http://getdp.info) libraries yourself, here are the options you will need.
The [PETSc](https://www.mcs.anl.gov/petsc) library should be compiled with the
following features:
--with-clanguage=cxx
--with-shared-libraries=1
--with-x=0
--with-scalar-type=complex
--with-mumps-serial=yes
--with-mpi=0
--with-mpiuni-fortran-binding=0
--download-mumps=yes
In other words, [PETSc](https://www.mcs.anl.gov/petsc) should be compiled
without [MPI](http://mpi-forum.org), with complex algebra and with a sequential
version of [MUMPS](http://mumps-solver.org). [GetDP](http://getdp.info) can then
be compiled with the following options:
ENABLE_BLAS_LAPACK: ON
ENABLE_BUILD_DYNAMIC: ON
ENABLE_BUILD_SHARED: ON
ENABLE_KERNEL: ON
ENABLE_PETSC: ON
ENABLE_WRAP_PYTHON: ON
If everything went fine, you should end up with the [GetDP](http://getdp.info)
Python module:
getdp.py
Finally, don't forget to update your `PATH`, `PYTHONPATH` and `LD_LIBRARY_PATH`.
Enjoy!
Reference paper
---------------
[W.-J., Beyn, "An integral method for solving nonlinear eigenvalue problems",
Linear Algebra and its Applications, 2012,436.](
https://doi.org/10.1016/j.laa.2011.03.030)
License
-------
Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
Institut fuer Theorie Elektromagnetischer Felder (TEMF),
Technische Universitaet Darmstadt.
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 2 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, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
cim-master/bin/000077500000000000000000000000001315622776500137055ustar00rootroot00000000000000cim-master/bin/args.py000066400000000000000000000062761315622776500152260ustar00rootroot00000000000000"""
cim.py, a non-linear eigenvalue solver.
Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
Institut fuer Theorie Elektromagnetischer Felder (TEMF),
Technische Universitaet Darmstadt.
See the LICENSE.txt and README.md for more license and copyright information.
"""
import argparse
def parse():
"""Parses python script arguments
Returns an argparse
"""
# Initialise
parser = argparse.ArgumentParser(prog="cim.py",
description=desc,
formatter_class=MyFormatter,
epilog="authors: " +
"N. Marsic, F. Wolf, " +
"S. Schoeps and H. De Gersem")
# Position argumets
parser.add_argument("pro", type=str,
help="GetDP .pro file")
parser.add_argument("mesh", type=str,
help="Gmsh .msh file")
parser.add_argument("resolution", type=str,
help="resolution from .pro file")
parser.add_argument("origin", type=complex,
help="circular contour origin (complex)")
parser.add_argument("radius", type=float,
help="circular contour radius")
# Optional arguments
parser.add_argument("-nodes", type=int, default=100,
help="number of nodes for trapezoidal rule")
parser.add_argument("-maxIt", type=int, default=10,
help="maximum number of iterations")
parser.add_argument("-lStart", type=int, default=4,
help="initial size of col(A0)")
parser.add_argument("-lStep", type=int, default=3,
help="step size for increasing col(A0)")
parser.add_argument("-rankTol", type=float, default=1e-4,
help="relative tolerance for rank test")
parser.add_argument("-setnumber", type=str, default=[], action='append',
nargs=2, metavar=('NAME', 'VALUE'),
help="set constant number NAME=VALUE")
parser.add_argument("-quiet", action='store_true',
help="should I be quiet?")
# Done
return parser.parse_args()
## Import only parse (other functions are just helpers)
__all__ = ['parse']
## Helper
desc = ("a non-linear eigenvalue solver " +
"using the contour integral method proposed by W.-J. Beyn " +
"coupled with a GetDP formulation.")
class MyFormatter(argparse.HelpFormatter):
""" Custom Formatter"""
def __init__(self, prog):
super(MyFormatter, self).__init__(prog,
indent_increment=2,
max_help_position=30,
width=None)
def _get_help_string(self, action):
help = action.help
if '%(default)' not in action.help:
if action.default is not argparse.SUPPRESS:
defaulting_nargs = [argparse.OPTIONAL, argparse.ZERO_OR_MORE]
if action.option_strings or action.nargs in defaulting_nargs:
help += ' (default: %(default)s)'
return help
cim-master/bin/beyn.py000066400000000000000000000163351315622776500152240ustar00rootroot00000000000000"""
cim.py, a non-linear eigenvalue solver.
Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
Institut fuer Theorie Elektromagnetischer Felder (TEMF),
Technische Universitaet Darmstadt.
See the LICENSE.txt and README.md for more license and copyright information.
"""
import solver as Solver
import sys
import numpy as np
import numpy.matlib
def simple(operator, origin, radius,
nNode=100, maxIt=10, lStart=4, lStep=3, rankTol=1e-4, verbose=True):
"""Solves an eigenvalue problem using Beyn's algorithm (simple version)
Keyword arguments:
operator -- the solver defining the operator to use
origin -- the origin (in the complex plane) of the above circular contour
radius -- the radius of the circular contour used to search the eigenvalues
nNode -- the number of nodes for the trapezoidal integration rule (optional)
maxIt -- the maximal number of iteration for constructing A0 (optional)
lStart -- the number of columns used for A0 when algorithm starts (optional)
lStep -- the step used for increasing the number of columns of A0 (optional)
rankTol -- the tolerance for the rank test (optional)
verbose -- should I be verbose? (optional)
Returns the computed eigenvalues, eigenvectors
and the associate *absolute* residual norms
"""
# Display the parameter used
if(verbose): display(nNode, maxIt, lStart, lStep, rankTol, origin, radius)
# Initialise A0 search
hasK = False
it = 0
m = operator.size()
l = lStart
k = -1
# Search A0
if(verbose): print "Searching A0 and A1 (" + str(m) + " unknowns)..."
while(not hasK and it != maxIt):
if(verbose): print " # Iteration " + str(it+1) + ": ",
if(verbose): sys.stdout.flush()
vHat = randomMatrix(m, l) # Take a random VHat
A0, A1 = integrate(operator, vHat, # Compute A_0 and A_1
nNode, radius, origin)
k, sV = rank(A0, rankTol) # Rank
if(verbose): print format(sV[0], ".2e") + " |",
if(verbose): print format(sV[-1], ".2e"),
if(verbose): print "[" + format(k, "d") + "]"
if(k == 0): # Rank is zero?
if(verbose): print " # Rank test is zero: keep A0 as is!"
hasK = True # -> Keep A0 (and warn)
elif(k == m): # Maximum rank reached?
if(verbose): print " # Maximum rank reached: keep A0 as is!"
hasK = True # -> Keep A0 (and warn)
elif(k < l): # Found a null SV?
hasK = True # -> We have A0
else: # Matrix is full rank?
l = l + lStep # -> Increase A0 size
it += 1 # Keep on searching A0
# Check if maxIt was reached
if(it == maxIt):
l = l - lStep
if(verbose): print "# Last iteration over: keep A0 as is!"
# Display singular values
if(verbose): print "Singular values: " + ", ".join("%.2e" % i for i in sV)
# Compute V, S and Wh
# NB: For SVD(A) = V*S*Wh, numpy computes {v, s, w}, such that:
# v = V; diag(s) = S and w = Wh
if(verbose): print "Constructing linear EVP..."
V, S, Wh = np.linalg.svd(A0, full_matrices=False, compute_uv=1)
# Extract V0, W0 and S0Inv
V0 = V[np.reshape(range(0, m), (-1, 1)), range(0, k)]
W0 = Wh.H[np.reshape(range(0, l), (-1, 1)), range(0, k)]
S0Inv = np.matrix(np.diag(1/S[0:k]))
# Compute B
B = V0.H * A1 * W0 * S0Inv
# Eigenvalues & eigenvectors
if(verbose): print "Solving linear EVP..."
myLambda, QHat = numpy.linalg.eig(B)
Q = V0 * QHat;
# Absolute residual norm
nFound = myLambda.shape[0]
norm = np.empty((nFound))
for i in range(nFound):
operator.apply(Q[:, i], myLambda[i]) # Apply eigenpair
r = operator.solution(0) # Get residual
norm[i] = np.linalg.norm(r) # Save the norm
# Done
if(verbose): print "Done!"
return myLambda, Q, norm
## Import only simple (other functions are just helpers)
__all__ = ["simple"]
## Helper functions
def rank(A, tol):
"""Returns the rank of a given matrix and its singular values
Keywords arguments:
A -- the matrix to use
tol -- the relative tolerance used to compute the rank
"""
S = np.linalg.svd(A, full_matrices=False, compute_uv=0)
nSV = len(S)
k = 0
for i in range(nSV):
if(S[i]/S[0] > tol): k = k + 1
return k, S.tolist()
def integrate(operator, vHat, nNode, radius, origin):
"""Computes the two first countour integrals of Beyn's method (A_0 and A_1)
over a circular contour.
Keyword arguments:
operator -- the solver defining the operator to use
vHat -- the RHS matrix defing Beyn's integrals
nNode -- the number of nodes used to discretise the circular contour
radius -- the radius of the circular contour
origin -- the origin (in the complex plane) of the circular contour
"""
# Initialise integrals
A0 = np.matlib.zeros(vHat.shape, dtype=complex)
A1 = np.matlib.zeros(vHat.shape, dtype=complex)
# Loop over integation points and integrate
for i in range(nNode):
t = 2 * np.pi * i / nNode;
phi = origin + radius * np.exp(1j * t)
tmp = multiSolve(operator, vHat, phi)
A0 += tmp * np.power(phi, 0) * np.exp(1j * t)
A1 += tmp * np.power(phi, 1) * np.exp(1j * t)
# Final scale
A0 *= radius/nNode
A1 *= radius/nNode
# Done
return A0, A1
def multiSolve(solver, B, w):
"""Solves for multiple RHS
Keyword arguments:
solver -- the solver to use
B -- the matrix of RHS
w -- the complex frequency to use
"""
# Solve
solver.solve(B, w)
# Get solutions
nSolve = B.shape[1]
size = solver.size()
X = np.matlib.empty((size, nSolve), dtype=complex)
for i in range(nSolve):
X[:, i] = solver.solution(i)
# Done
return X
def randomMatrix(n, m):
"""Returns a random complex matrix of the given size (n x m)"""
return np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j
def display(nNode, maxIt, lStart, lStep, rankTol, origin, radius):
print "Beyn's contour integral method (simple)"
print "---------------------------------------"
print " # Nodes used for the trapezoidal rule:" + " " + str(nNode)
print " # Maximum number of iterations: " + " " + str(maxIt)
print " # Initial size of col(A0): " + " " + str(lStart)
print " # Step size for increasing col(A0): " + " " + str(lStep)
print " # Rank test relative tolerance: ",
print format(rankTol, ".2e")
print "---------------------------------------"
print " # Cirular path origin: ",
print "(" + format(origin.real, "+.2e") + ")",
print "+",
print "(" + format(origin.imag, "+.2e") + ")j"
print " # Cirular path radius: ",
print format(radius, "+.2e")
print "---------------------------------------"
cim-master/bin/cim.py000077500000000000000000000034731315622776500150410ustar00rootroot00000000000000#!/usr/bin/env python2
"""
cim.py, a non-linear eigenvalue solver.
Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
Institut fuer Theorie Elektromagnetischer Felder (TEMF),
Technische Universitaet Darmstadt.
See the LICENSE.txt and README.md for more license and copyright information.
"""
import solver as Solver
import beyn as Beyn
import numpy as np
import args as args
import time as timer
# Start timer
tStart = timer.time()
# Parse arguments
arg = args.parse()
# Prepare variables of -setnumber for GetDP
setnumber = []
for i in arg.setnumber:
setnumber.append("-setnumber")
for j in i:
setnumber.append(j)
# Initialise GetDP
operator = Solver.GetDPWave(arg.pro, arg.mesh, arg.resolution, setnumber)
# Compute
l, V, r = Beyn.simple(operator, arg.origin, arg.radius,
arg.nodes, arg.maxIt, arg.lStart, arg.lStep, arg.rankTol,
not arg.quiet)
# Display the computed eigenvalues
if(not arg.quiet):
print
print "Eigenvalues:"
print "----------- "
for i in range(l.shape[0]):
if(not arg.quiet): print " # " + str(i) + ": ",
print "(" + format(np.real(l[i]).tolist(), '+.12e') + ")",
print "+",
print "(" + format(np.imag(l[i]).tolist(), '+.12e') + ")*j",
print "| " + format(r[i], 'e'),
print "| " + format(r[i]/np.abs(l[i]), 'e')
if(not arg.quiet):
print "----------- "
print "Found: " + str(l.shape[0])
# Generate GetDP views for eigenvectors
if(not arg.quiet):
print
print "Generating eigenvector views..."
for i in range(l.shape[0]):
operator.view(V[:, i], "cim" + str(i) + ".pos")
if(not arg.quiet): print " # View: " + str(i) + " done!"
# Stop timer
tStop = timer.time()
# Done
if(not arg.quiet):
print
print "Elapsed time: " + str(round(tStop - tStart)) + "s"
print "Bye!"
cim-master/bin/solver.py000066400000000000000000000114111315622776500155670ustar00rootroot00000000000000"""
cim.py, a non-linear eigenvalue solver.
Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
Institut fuer Theorie Elektromagnetischer Felder (TEMF),
Technische Universitaet Darmstadt.
See the LICENSE.txt and README.md for more license and copyright information.
"""
import getdp as GetDP
import numpy as np
class GetDPWave:
"""A GetDP solver using complex arithmetic for time-harmonic wave problems
Only one instance of this class makes sens
The .pro file should use the following variables:
angularFreqRe -- the real part of the angular frequency to use
angularFreqIm -- the imaginary part of the angular frequency to use
nRHS -- the number of right hand side that should be considered (default: 1)
b_i / b~{i} -- the ith right hand side (RHS) to use (first index is 0)
x_i / x~{i} -- the ith solution vector computed by GetDP (first index is 0)
doInit -- should some initialization be done (default: 0)?
doSolve -- should Ax_i = b_i be solved for all i (default: 0)?
doPostpro -- should only a view be created for x (default: 0)?
doApply -- should only the application of x be computed (default: 0)?
fileName -- post-processing file name
"""
def __init__(self, pro, mesh, resolution, optional=[]):
"""Instanciates a new SolverGetDP with a full '-solve'
Keyword arguments:
pro -- the .pro file to use
mesh -- the .msh file to use
resolution -- the resolution (from the .pro file) to use
optional -- optional arguments for GetDP (default value = [])
Optional argument has the following structure:
["GetDP option", "value 1", "value 2", ..., "GetDP option", ...]
"""
# Save
self.__pro = pro
self.__mesh = mesh
self.__resolution = resolution
self.__optional = optional
# Generate DoFs and assemble a first system
GetDP.GetDPSetNumber("doInit", 1);
GetDP.GetDP(["getdp", self.__pro,
"-msh", self.__mesh,
"-solve", self.__resolution,
"-v", "2"] + self.__optional)
def solution(self, i):
"""Returns the solution for the ith RHS (first index is 0)"""
return self.__toNumpy(GetDP.GetDPGetNumber("x_" + str(i)))
def size(self):
"""Returns the number of degrees of freedom"""
return self.solution(0).shape[0]
def apply(self, x, w):
"""Applies x to the operator with a pulsation of w
This method updates self.solution()
"""
GetDP.GetDPSetNumber("angularFreqRe", np.real(w).tolist())
GetDP.GetDPSetNumber("angularFreqIm", np.imag(w).tolist())
GetDP.GetDPSetNumber("x_0", self.__toGetDP(x))
GetDP.GetDPSetNumber("doApply", 1)
GetDP.GetDP(["getdp", self.__pro,
"-msh", self.__mesh,
"-cal"] + self.__optional)
def solve(self, b, w):
"""Solves with b as RHS and w as complex angular frequency
b is a matrix and each column is a different RHS
"""
# Number of RHS
nRHS = b.shape[1]
# Set variables
GetDP.GetDPSetNumber("nRHS", nRHS)
GetDP.GetDPSetNumber("angularFreqRe", np.real(w).tolist())
GetDP.GetDPSetNumber("angularFreqIm", np.imag(w).tolist())
# Set RHS
for i in range(nRHS):
GetDP.GetDPSetNumber("b_" + str(i), self.__toGetDP(b[:, i]))
# Solve
GetDP.GetDPSetNumber("doSolve", 1)
GetDP.GetDP(["getdp", self.__pro,
"-msh", self.__mesh,
"-cal"] + self.__optional)
def view(self, x, fileName):
"""Generates a post-processing view
Keyword arguments:
x -- the solution vector to use
fileName -- the post-precessing file name
This method generates a linear system
"""
GetDP.GetDPSetNumber("x_0", self.__toGetDP(x))
GetDP.GetDPSetNumber("doPostpro", 1)
GetDP.GetDPSetString("fileName", fileName)
GetDP.GetDP(["getdp", self.__pro,
"-msh", self.__mesh,
"-cal"] + self.__optional)
@staticmethod
def __toNumpy(vGetDP):
"""Takes a GetDP list and returns a numpy array"""
size = vGetDP.size() / 2
vNumpy = np.empty((size, 1), dtype=complex)
for i in range(size):
vNumpy[i] = complex(vGetDP[i*2], vGetDP[i*2 + 1])
return vNumpy
@staticmethod
def __toGetDP(vNumpy):
"""Takes a numpy array and returns a GetDP list"""
size = vNumpy.shape[0]
vGetDP = list()
for i in range(size):
vGetDP.append(float(np.real(vNumpy[i])))
vGetDP.append(float(np.imag(vNumpy[i])))
return vGetDP
cim-master/dependencies.sh000077500000000000000000000042361315622776500161270ustar00rootroot00000000000000#!/bin/sh
## Dependencies folder ##
#########################
echo "Creating dependencies folder: dep/"
mkdir dep
cd dep
## BLAS ##
##########
echo "Downloading OpenBLAS"
wget http://github.com/xianyi/OpenBLAS/archive/v0.2.19.tar.gz
echo "Extracting OpenBLAS"
tar xf v0.2.19.tar.gz
echo "Entering OpenBLAS"
cd OpenBLAS-0.2.19
BLAS=$(pwd)
echo "Compiling OpenBLAS"
make
echo "Leaving OpenBLAS"
cd ..
## PETSc ##
###########
echo "Downloading PETSc"
wget http://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-3.7.5.tar.gz
echo "Extracting PETSc"
tar xf petsc-3.7.5.tar.gz
echo "Entering PETSc"
cd petsc-3.7.5
PETSC_DIR=$(pwd)
PETSC_ARCH=linux-gnu-cxx-seq
PETSC=$PETSC_DIR"/"$PETSC_ARCH
echo "Configuring PETSc"
./configure --with-clanguage=cxx --with-shared-libraries=1 --with-x=0 \
--with-scalar-type=complex --with-mumps-serial=yes --with-mpi=0 \
--with-mpiuni-fortran-binding=0 --download-mumps=yes \
--with-blas-lapack-lib=$BLAS/libopenblas.so
echo "Compiling PETSc"
make all
echo "Checking PETSc"
make test
echo "Leaving PETSc"
cd ..
## GetDP ##
###########
echo "Downloading GetDP"
wget http://getdp.info/src/getdp-2.11.0-source.tgz
echo "Extracting GetDP"
tar xf getdp-2.11.0-source.tgz
echo "Entering GetDP"
cd getdp-2.11.0-source
echo "Configuring GetDP"
mkdir build
cd build
cmake -DDEFAULT=0 -DENABLE_BUILD_DYNAMIC=1 -DENABLE_BUILD_SHARED=1 \
-DENABLE_KERNEL=1 -DENABLE_PETSC=1 -DENABLE_WRAP_PYTHON=1 \
-DENABLE_BLAS_LAPACK=1 -DBLAS_LAPACK_LIBRARIES=$BLAS/libopenblas.so \
-DPython_ADDITIONAL_VERSIONS="2.7" ..
echo "Compiling GetDP"
GETDP=$(pwd)
make
echo "Leaving GetDP"
cd ../../..
## PATH ##
##########
echo " "
echo "+++"
echo "PETSc and GetDP should be now compiled with the required features."
echo "In order to use cim.py, a few additional steps are required."
echo "We still need to add the libraries and the Python module in your path."
echo "Please add the following in your ~/.bashrc file (or equivalent):"
echo "export PATH=\$PATH:"$(pwd)"/bin"
echo "export PYTHONPATH=\$PYTHONPATH:"$GETDP
echo "export LD_LIBRARY_PATH=\$LD_LIBRARY_PATH:"$PETSC"/lib:"$GETDP":"$BLAS
echo "+++"
## Enjoy ! ##
#############
cim-master/example/000077500000000000000000000000001315622776500145705ustar00rootroot00000000000000cim-master/example/maxwell_linear/000077500000000000000000000000001315622776500175735ustar00rootroot00000000000000cim-master/example/maxwell_linear/Makefile000066400000000000000000000005601315622776500212340ustar00rootroot00000000000000BIN=../../bin
test:
@echo "make cim: run the countour integral solver"
@echo "make ref: compute the first eigenvalue directly with GetDP"
ref: init
getdp ref.pro -solve Eig -pos Eig -msh square.msh
cim: init
python ${BIN}/cim.py maxwell.pro square.msh Maxwell 9e8 1e8 -lStart 1 -lStep 1
init:
gmsh square.geo -2
clean:
rm -f *.pos
rm -f *.pre
rm -f *.msh
cim-master/example/maxwell_linear/README.md000066400000000000000000000013361315622776500210550ustar00rootroot00000000000000This is a simple linear eigenvalue problem example. It consists in a rectangular
electromagnetic cavity. Simply run:
make cim
to test cim.py.
The [Makefile](Makefile) will first mesh the geometry [square.geo](square.geo)
by calling [Gmsh](http://gmsh.info).
Afterwards, cim.py is called. The [GetDP](http://getdp.info) formulation is
located in [maxwell.pro](maxwell.pro). In order to check the solution,
simply run:
make ref
This will solve the linear eigenvalue problem with a classical algorithm. This
resolution is implemented in [ref.pro](ref.pro). If you use [GetDP](
http://getdp.info) with [SLEPc](http://slepc.upv.es), the default algorithm
should be [Krylov-Schur](https://dx.doi.org/10.1137/S0895479800371529).
cim-master/example/maxwell_linear/maxwell.pro000066400000000000000000000062311315622776500217700ustar00rootroot00000000000000Group{
Omega = Region[1]; // Domain
Bndry = Region[2]; // Boundary
}
Function{
// Physical data //
DefineConstant[angularFreqRe = 5, // [rad/s]
angularFreqIm = 0]; // [rad/s]
C0 = 299792458; // [m/s]
aFC[] = Complex[angularFreqRe,
angularFreqIm]; // [rad/m]
k[] = aFC[]/C0; // [rad/m]
Printf["Angular frequency set to: %g + i*%g [rad/s]",
angularFreqRe, angularFreqIm];
// Algebraic data //
DefineConstant[nRHS = 1]; // Number of RHS for this run
For i In {0:nRHS-1}
DefineConstant[x~{i}() = {}, // Solution
b~{i}() = {}]; // Right hand side
EndFor
// Control data //
DefineConstant[doInit = 0, // Should I initialize some stuff?
doSolve = 0, // Should I solve Ax = b?
doPostpro = 0, // Should I only create a view for x()?
doApply = 0, // Should I only apply x(): x <- Ax?
fileName = "eig.pos"]; // Postpro file name
}
Jacobian{
{ Name JVol;
Case{
{ Region All; Jacobian Vol; }
}
}
}
Integration{
{ Name IP2;
Case{
{ Type Gauss;
Case{
{ GeoElement Line; NumberOfPoints 2; }
{ GeoElement Triangle; NumberOfPoints 3; }
{ GeoElement Tetrahedron; NumberOfPoints 4; }
}
}
}
}
}
Constraint{
{ Name Dirichlet; Type Assign;
Case{
{ Region Bndry; Value 0; }
}
}
}
FunctionSpace{
{ Name HCurl; Type Form1;
BasisFunction{ // Nedelec
{ Name se; NameOfCoef ee; Function BF_Edge;
Support Omega; Entity EdgesOf[All]; }
}
Constraint{
{ NameOfCoef ee; EntityType EdgesOf; NameOfConstraint Dirichlet; }
}
}
}
Formulation{
{ Name Maxwell; Type FemEquation;
Quantity{
{ Name e; Type Local; NameOfSpace HCurl; }
}
Equation{
Galerkin{ [ Dof{d e} , {d e}];
In Omega; Integration IP2; Jacobian JVol; }
Galerkin{ [-k[]^2 * Dof{ e} , { e}];
In Omega; Integration IP2; Jacobian JVol; }
}
}
}
Resolution{
{ Name Maxwell;
System{
{ Name A; NameOfFormulation Maxwell; Type Complex; }
}
Operation{
Generate[A];
If(doInit)
CopySolution[A, x~{0}()];
EndIf
If(doSolve)
// Full solve for first RHS
CopyRightHandSide[b~{0}(), A];
Solve[A];
CopySolution[A, x~{0}()];
// SolveAgain for other RHSs
For i In {1:nRHS-1}
CopyRightHandSide[b~{i}(), A];
SolveAgain[A];
CopySolution[A, x~{i}()];
EndFor
EndIf
If(doApply)
CopySolution[x~{0}(), A];
Apply[A];
CopySolution[A, x~{0}()];
EndIf
If(doPostpro)
CopySolution[x~{0}(), A];
PostOperation[Maxwell];
EndIf
}
}
}
PostProcessing{
{ Name Maxwell; NameOfFormulation Maxwell;
Quantity{
{ Name e; Value{ Local{ [{e}]; In Omega; Jacobian JVol; } } }
}
}
}
PostOperation{
{ Name Maxwell; NameOfPostProcessing Maxwell;
Operation{
Print[e, OnElementsOf Omega, File fileName];
}
}
}
cim-master/example/maxwell_linear/ref.pro000066400000000000000000000034461315622776500211000ustar00rootroot00000000000000Group{
Omega = Region[1]; // Domain
Bndry = Region[2]; // Boundary
}
Function{
// Physical data //
C0 = 299792458; // [m/s]
// Eigenvalue solver data //
DefineConstant[nEig = 10, // [-]
target = 8e17]; // [(rad/s)^2]
}
Jacobian{
{ Name JVol;
Case{
{ Region All; Jacobian Vol; }
}
}
}
Integration{
{ Name IP2;
Case{
{ Type Gauss;
Case{
{ GeoElement Line; NumberOfPoints 2; }
{ GeoElement Triangle; NumberOfPoints 3; }
{ GeoElement Tetrahedron; NumberOfPoints 4; }
}
}
}
}
}
Constraint{
{ Name Dirichlet; Type Assign;
Case{
{ Region Bndry; Value 0; }
}
}
}
FunctionSpace{
{ Name HCurl; Type Form1;
BasisFunction{ // Nedelec
{ Name se; NameOfCoef ee; Function BF_Edge;
Support Omega; Entity EdgesOf[All]; }
}
Constraint{
{ NameOfCoef ee; EntityType EdgesOf; NameOfConstraint Dirichlet; }
}
}
}
Formulation{
{ Name Eig; Type FemEquation;
Quantity{
{ Name e; Type Local; NameOfSpace HCurl; }
}
Equation{
Galerkin{ [ Dof{d e} , {d e}];
In Omega; Integration IP2; Jacobian JVol; }
Galerkin{ DtDtDof[1/C0^2 * Dof{ e} , { e}];
In Omega; Integration IP2; Jacobian JVol; }
}
}
}
Resolution{
{ Name Eig;
System{
{ Name A; NameOfFormulation Eig; Type Complex; }
}
Operation{
GenerateSeparate[A];
EigenSolve[A, nEig, target, 0];
}
}
}
PostProcessing{
{ Name Eig; NameOfFormulation Eig;
Quantity{
{ Name e; Value{ Local{ [{e}]; In Omega; Jacobian JVol; } } }
}
}
}
PostOperation{
{ Name Eig; NameOfPostProcessing Eig;
Operation{
Print[e, OnElementsOf Omega, File "eig.pos", EigenvalueLegend];
}
}
}
cim-master/example/maxwell_linear/square.geo000066400000000000000000000010531315622776500215660ustar00rootroot00000000000000DefineConstant[l = {01, Name "Input/00Geometry/00Length [m]"},
msh = {10, Name "Input/01Mesh/00Density [point per length]"}];
Point(1) = {+l/2, -l/2, 0, 1};
Point(2) = {+l/2, +l/2, 0, 1};
Point(3) = {-l/2, +l/2, 0, 1};
Point(4) = {-l/2, -l/2, 0, 1};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Transfinite Line {1, 2, 3, 4} = (msh * l) Using Progression 1;
Transfinite Surface {1};
Physical Surface(1) = {1};
Physical Line(2) = {1, 2, 3, 4};
cim-master/example/maxwell_sibc/000077500000000000000000000000001315622776500172415ustar00rootroot00000000000000cim-master/example/maxwell_sibc/Makefile000066400000000000000000000003561315622776500207050ustar00rootroot00000000000000BIN=../../bin
sphere: init cim
init:
gmsh sphere.geo -3
cim:
python ${BIN}/cim.py sphere.pro sphere.msh Solve 8.22543e9+2.46361e1j 1e8 -nodes 10 -lStart 4 -setnumber isSC 0
clean:
rm -f *.pos
rm -f *.msh
rm -f *.pre
rm -f *.res
cim-master/example/maxwell_sibc/README.md000066400000000000000000000021271315622776500205220ustar00rootroot00000000000000This is a non-linear eigenvalue problem example. It consists in a spherical
electromagnetic cavity. The wall conductivity is modelled by a Leontovich
surface impedance boundary condition (SIBC). Simply run:
make
to test cim.py.
The [Makefile](Makefile) will first mesh the geometry [sphere.geo](sphere.geo)
by calling [Gmsh](http://gmsh.info). Afterwards, cim.py is called. The [GetDP](
http://getdp.info) formulation is located in [sphere.pro](sphere.pro).
The [GetDP](http://getdp.info) codes
related to cim.py are located in [cimParameters.pro](cimParameters.pro) and
[cimResolution.pro](cimResolution.pro).
For the default parameters:
radius = 100 mm
conductivity = 1e15 S/m
the analytical resonance angular frequency for the fundamental mode should be
8.22543e9+2.46361e1j
where j is the imaginary unit. This analytical result comes from the reference:
[S. Papantonis and S. Lucyszyn, "Lossy spherical cavity resonators for
stress-testing arbitrary 3D eigenmode solvers," Progress In Electromagnetics
Research, vol. 151, pp. 151-167, 2015.](https://dx.doi.org/10.2528/PIER15031702)
cim-master/example/maxwell_sibc/cimParameters.pro000066400000000000000000000020721315622776500225600ustar00rootroot00000000000000/////////////////////////////////////////////////////////////////
// Parameters for contour integral method. //
// The master file should inclue "cimResolution.pro" later on. //
/////////////////////////////////////////////////////////////////
Function{
// Physical data //
DefineConstant[angularFreqRe = 1, // [rad/s]
angularFreqIm = 0]; // [rad/s]
Puls[] = Complex[angularFreqRe,
angularFreqIm]; // [rad/m]
// Algebraic data //
DefineConstant[nRHS = 1]; // Number of RHS for this run
For i In {0:nRHS-1}
DefineConstant[x~{i}() = {}, // Solution
b~{i}() = {}]; // Right hand side
EndFor
// Control data //
DefineConstant[doInit = 0, // Should I initialize some stuff?
doSolve = 0, // Should I solve Ax = b?
doPostpro = 0, // Should I only create a view for x()?
doApply = 0, // Should I only apply x(): x <- Ax?
fileName = "eig.pos"]; // Postpro file name
}
cim-master/example/maxwell_sibc/cimResolution.pro000066400000000000000000000021411315622776500226150ustar00rootroot00000000000000/////////////////////////////////////////////////////////////////////
// Resolution for contour integral method. //
// Warning, the master file should include "cimParameters.pro", //
// and should also define "FormulationCIM" the formulation to use. //
/////////////////////////////////////////////////////////////////////
Resolution{
{ Name Solve;
System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
Operation{
Generate[A];
If(doInit)
CopySolution[A, x~{0}()];
EndIf
If(doSolve)
// Full solve for first RHS
CopyRightHandSide[b~{0}(), A];
Solve[A];
CopySolution[A, x~{0}()];
// SolveAgain for other RHSs
For i In {1:nRHS-1}
CopyRightHandSide[b~{i}(), A];
SolveAgain[A];
CopySolution[A, x~{i}()];
EndFor
EndIf
If(doApply)
CopySolution[x~{0}(), A];
Apply[A];
CopySolution[A, x~{0}()];
EndIf
If(doPostpro)
CopySolution[x~{0}(), A];
PostOperation[Post];
EndIf
}
}
}
cim-master/example/maxwell_sibc/sphere.dat000066400000000000000000000014111315622776500212160ustar00rootroot00000000000000// Geometry //
DefineConstant[R = { 100e-3, Name "Input/00Geometry/00Radius [m]" }];
// Mesh //
DefineConstant[md = { 10, Name "Input/01Mesh/00Density [Radius^-1]" }];
cl = R/md;
// Material //
DefineConstant[isSC = { 0,
Name "Input/02Material/00Is superconductor?",
GmshOption "Reset",
Choices {0 = "No", 1 = "Yes"} }];
If(isSC == 0)
DefineConstant[Sigma = { 1e15,
Name "Input/02Material/01Conductivity [S m^-1]" }];
EndIf
If(isSC == 1)
DefineConstant[Sigma = { 1e9,
Name "Input/02Material/01Conductivity [S m^-1]" }];
DefineConstant[Lambda = { 1e-8,
Name "Input/02Material/02London depth [m]" }];
EndIf
cim-master/example/maxwell_sibc/sphere.geo000066400000000000000000000021121315622776500212170ustar00rootroot00000000000000Include "sphere.dat";
Point(0) = { 0, 0, 0, cl};
Point(1) = {+R, 0, 0, cl};
Point(2) = { 0, +R, 0, cl};
Point(3) = {-R, 0, 0, cl};
Point(4) = { 0, -R, 0, cl};
Point(5) = { 0, 0, +R, cl};
Point(6) = { 0, 0, -R, cl};
Circle(01) = {1, 0, 2};
Circle(02) = {2, 0, 3};
Circle(03) = {3, 0, 4};
Circle(04) = {4, 0, 1};
Circle(05) = {5, 0, 1};
Circle(06) = {1, 0, 6};
Circle(07) = {6, 0, 3};
Circle(08) = {3, 0, 5};
Circle(09) = {2, 0, 6};
Circle(10) = {6, 0, 4};
Circle(11) = {4, 0, 5};
Circle(12) = {5, 0, 2};
Line Loop(1) = {1, -12, 5};
Line Loop(2) = {1, 9, -6};
Line Loop(3) = {9, 7, -2};
Line Loop(4) = {2, 8, 12};
Line Loop(5) = {11, 5, -4};
Line Loop(6) = {4, 6, 10};
Line Loop(7) = {10, -3, -7};
Line Loop(8) = {3, 11, -8};
Ruled Surface(1) = {1};
Ruled Surface(2) = {2};
Ruled Surface(3) = {3};
Ruled Surface(4) = {4};
Ruled Surface(5) = {5};
Ruled Surface(6) = {6};
Ruled Surface(7) = {7};
Ruled Surface(8) = {8};
Surface Loop(1) = {3, 2, 1, 4, 8, 7, 6, 5};
Volume(1) = {1};
Physical Volume(1) = {1};
Physical Surface(2) = {1, 2, 3, 4, 5, 6, 7, 8};
Physical Point(3) = {0};
cim-master/example/maxwell_sibc/sphere.pro000066400000000000000000000041251315622776500212530ustar00rootroot00000000000000Include "sphere.dat";
Group{
Omega = Region[1];
Bndry = Region[2];
}
Include "./cimParameters.pro"; // Functions for cim.py
Function{
// Imaginary Unit //
J[] = Complex[0, 1];
// Physical Constants //
C0 = 299792458; // [m/s]
Mu0 = 4*Pi*1e-7; // [H/m]
Eps0 = 1/(Mu0 * C0^2); // [F/m]
// Leontovich SIBC (H-Formulation) //
If(isSC == 0)
SL[] = Sqrt[(J[] * Puls[] * Mu0) / Sigma]; // H-Formulation
EndIf
If(isSC == 1)
SL[] = 0.5 * Puls[]^2 * Mu0^2 * Lambda^3 * Sigma +
J[] * Puls[] * Mu0 * Lambda; // H-Formulation
EndIf
SIbc[] = -1.0 / SL[]; // E-Formulation
}
Integration{
{ Name I1;
Case{
{ Type Gauss;
Case{
{ GeoElement Tetrahedron; NumberOfPoints 4; } // O1 * O1
{ GeoElement Triangle; NumberOfPoints 3; } // O1 * O1
{ GeoElement Line; NumberOfPoints 2; } // O1 * O1
}
}
}
}
}
Jacobian{
{ Name Jac;
Case{
{ Region Omega; Jacobian Vol; }
{ Region Bndry; Jacobian Sur; }
}
}
}
FunctionSpace{
{ Name HCurl; Type Form1;
BasisFunction{
{ Name sx; NameOfCoef ux; Function BF_Edge;
Support Region[{Omega, Bndry}]; Entity EdgesOf[All]; }
}
}
}
Formulation{
{ Name FormulationCIM; Type FemEquation;
Quantity{
{ Name e; Type Local; NameOfSpace HCurl; }
}
Equation{
// Maxwell
Galerkin{ [ 1/Mu0 * Dof{d e}, {d e}];
In Omega; Jacobian Jac; Integration I1; }
Galerkin{ [ -Puls[]^2 * Eps0 * Dof{ e}, { e}];
In Omega; Jacobian Jac; Integration I1; }
// IBC
Galerkin{ [-J[] * Puls[] * SIbc[] * Dof{ e}, { e}];
In Bndry; Jacobian Jac; Integration I1; }
}
}
}
Include "./cimResolution.pro"; // Resolution for cim.py
PostProcessing{
{ Name Post; NameOfFormulation FormulationCIM;
Quantity{
{ Name E; Value{ Term{ [{e}]; In Omega; Jacobian Jac; } } }
}
}
}
PostOperation{
{ Name Post; NameOfPostProcessing Post;
Operation{
Print[E, OnElementsOf Omega, File fileName];
}
}
}